A geometrically motivated coordinate system for exploring spacetime dynamics in 
numerical-relativity simulations using a quasi-Kinnersley tetrad 
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We investigate the suitability and properties of a quasi-Kinnersley tetrad and a geometrically 
motivated coordinate system as tools for quantifying both strong-field and wave-zone effects in nu- 
merical relativity (NR) simulations. We fix two of the coordinate degrees of freedom of the metric, 
namely the radial and latitudinal coordinates, using the Coulomb potential associated with the 
quasi-Kinnersley transverse frame. These coordinates are invariants of the spacetime and can be 
used to unambiguously fix the outstanding spin-boost freedom associated with the quasi-Kinnersley 
frame (and thus can be used to choose a preferred quasi-Kinnersley tetrad). In the limit of small per- 
turbations about a Kerr spacetime, these geometrically motivated coordinates and quasi-Kinnersley 
tetrad reduce to Boyer-Lindquist coordinates and the Kinnersley tetrad, irrespective of the simula- 
tion gauge choice. We explore the properties of this construction both analytically and numerically, 
and we gain insights regarding the propagation of radiation described by a super-Poynting vector, 
further motivating the use of this construction in NR simulations. We also quantify in detail the 
peeling properties of the chosen tetrad and gauge. We argue that these choices are particularly 
well suited for a rapidly converging wave-extraction algorithm as the extraction location approaches 
infinity, and we explore numerically the extent to which this property remains applicable on the 
interior of a computational domain. Using a number of additional tests, we verify numerically that 
the prescription behaves as required in the appropriate limits regardless of simulation gauge; these 
tests could also serve to benchmark other wave extraction methods. We explore the behavior of the 
geometrically motivated coordinate system in dynamical binary-black-hole NR mergers; while we 
obtain no unexpected results, we do find that these coordinates turn out to be useful for visualizing 
NR simulations (for example, for vividly illustrating effects such as the initial burst of spurious 
"junk" radiation passing through the computational domain). Finally, we carefully scrutinize the 
head-on collision of two black holes and, for example, the way in which the extracted waveform 
changes as it moves through the computational domain. 

PACS numbers: 04.25.D-,04.30.-w,04.25.dg 



I. INTRODUCTION 

Numerical relativity (NR) has made great strides in 
recent years and is now able to explore binary black 
hole, black hole - neutron star, and neutron star - neu- 
tron star mergers in a wide variety of configurations 
(see [1, 2] for recent reviews). Numerical simulations 
are crucial tools for calibrating and validating the large 
template bank of analytic, phenomenological waveforms 
that will be used to search for gravitational waves in 
data from detectors such as the Laser Interferometer 
Gravitational- Wave Observatory (LIGO) [3, 4], Virgo [5] 
and the Large-scale Cryogenic Gravitational-wave Tele- 
scope (LCGT) [6]. Numerical simulations also make it 
possible, for the first time, to explore fully dynamical 
spacetimes in the strong field region, such as the space- 
time of a compact-binary merger. 

An important attribute of any analysis performed on 
numerical simulations is the ability to extract informa- 
tion in a manner independent of the gauge in which one 
chooses to perform the simulation. In this paper, we sug- 
gest one such strategy: using a quasi-Kinnersley tetrad 
adapted to a choice of coordinates that are computed us- 
ing the curvature invariants of the spacetime. Most cal- 



culations presented in this paper are local, allowing our 
tetrad and choice of geometrically motivated coordinates 
(and all quantities derived from them) to be computed 
in real time during a numerical simulation (i.e. without 
post-processing). The proposed scheme is also applica- 
ble throughout the spacetime, allowing us to study phe- 
nomena in both the strongly curved and asymptotic flat 
regions with the same tools. 

In order to extract the 6 physical degrees of freedom of 
a general Lorentzian metric in four dimensions expressed 
in terms of a tetrad formulation, 10 degrees of freedom 
have to be specified [see e.g. [7]]. Of these 10 degrees of 
freedom, 6 are associated with the tetrad at a particular 
point on the spacetime manifold and 4 originate from 
the freedom to label that point (the choice of gauge). 
A common choice of tetrad and the one adopted here 
is the Newman-Penrosc (NP) null tetrad, which consists 
of two real null vectors denoted I and n as well as a 
complex conjugate pair of null vectors m and m. As we 
demonstrate explicitly in Sec. II, where we consider the 
mathematical details in greater depth, the tetrad choice 
is not unique. The freedom to orient and scale the tetrad 
is expressed by 6 parameters associated with a general 
Lorentz transformation between two different null tetrads 
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at a fixed point in the spacctime. 

In order to extract physically meaningful quantities 
and to compare results from different simulations and 
numerical codes, an explicit prescription for the tetrad 
choice must be made. Two geometrically motivated pre- 
scriptions for orientating the tetrad immediately suggest 
themselves: choosing i) the principle null frame (PNF) 
or ii) the transverse frame (TF) . (By "frame" , we mean 
a set of tetrads related by a Type III transformation 
[Sec. IIC].) The relationship between these two frames 
and their properties are discussed in greater detail in 
Sees. II and III; either one of these two choices imme- 
diately removes 4 of the 6 possible tetrad degrees of free- 
dom. The remaining 2 degrees of freedom in the tetrad 
choice are more subtle [see the discussion in Sec. IIIC]. 

The procedure we adopt in this paper is to choose 
a special transverse tetrad that becomes the Kinnersley 
tetrad [8] in Type-D spacetimes. The properties of these 
tetrads (known as quasi-Kinnersley tetrads, or QKTs) 
and their importance for NR have previously been ex- 
plored in detail [9-14] . Part of the motivation for choos- 
ing a QKT is implicit in Chandrasekhar's work on the 
gravitational perturbations of the Kerr black hole [15]: 
in this work, he showed that for a given perturbation of 
the Kerr metric (expressed in terms of the Weyl scalar 
5^4) it is always possible, working to linear order, to 
find a transverse tetrad and a gauge constructed from 
the Coulomb potential associated with this tetrad such 
that the Coulomb potential of the perturbed and back- 
ground metrics are the same. 

We revisit these ideas in Sec. Ill, where we investigate 
the properties of the quasi-Kinnersley tetrad choice. We 
concentrate on the implications of the intrinsic geomet- 
rical properties of the tetrad, rather than (as previous 
works have done) focusing on the tetrad's properties in 
a perturbative limit. For example, we explore the di- 
rections of energy flow using the super-Poynting vector, 
showing that the choice of a QKT naturally aligns the 
tetrad with the wave-fronts of passing radiation. This ob- 
servation suggests that, even in the strong field regime, 
the QKT is a natural, geometrically motivated tetrad 
choice for observing the flow of radiation and other space- 
time dynamics. 

After specializing to the transverse frame there exist 
two remaining degrees of tetrad freedom: the freedom 
of the spin-boost transformations. We fix this remaining 
tetrad freedom by relying on a straightforward extension 
of Chandrasekhar's work [15] to the strong field regime. 
We present the mathematical details in Sees. HID and 
HIE: specifically, we use the Coulomb potential ^2 on 
the QKT to introduce a pair of geometrically motivated 
radial and latitudinal coordinates. Note that ^2 on the 
transverse frame is spin-boost independent, that the re- 
sulting coordinates can be constructed from the curva- 
ture invariants I and J only, and that these coordinates 
reduce to the Boycr-Lindquist radial and latitudinal co- 
ordinates for Kerr spacetimes. We then use these ge- 
ometrically motivated coordinates to fix the spin-boost 



freedom by ensuring that the projection of the tetrad 
base vectors onto the gradients of the new coordinates 
obey the relations found for the Kinnersley tetrad in the 
Kerr limit. 

One application of our chosen QKT and geometrically 
motivated coordinates is gravitational-wave extraction. 
For isolated, gravitating systems, gravitational radiation 
is only strictly defined at future null infinity; this is a con- 
sequence of the so-called "peeling property" that governs 
the behavior of the Weyl curvature scalars as measured 
on an affinely parametrized out-going geodesic. With 
the goal of using our tetrad and gauge prescription as a 
possible wave extraction method, we work out the impli- 
cations that this peeling property has for the Weyl curva- 
ture scalar expressed on the QKT in Sec. IIIH. We high- 
light not only the falloff behavior of the QKT Ncwman- 
Penrose quantities but also the behavior of the geomet- 
rically constructed radial coordinate. We explore graph- 
ically some of the implications of the peeling property 
for the bunching of principle null directions and argue 
that the directions associated with QKT are the optimal 
out-going directions for ensuring rapid convergence of the 
computed radiation quantities to the correct asymptotic 
results. 

We implement our geometrically motivated coordi- 
nates and QKT numerically within the context of a pseu- 
dospectral NR code in Sec. IV, and we present a number 
of numerical simulations demonstrating the behavior of 
our coordinates and QKT in Sec. V. First, we carry out, 
for both non-radiative and radiative spacetimes, a few 
checks to verify that our scheme works correctly regard- 
less of the choice of gauge in the simulation [Sees. VA 
and VB, respectively]. We confirm that we obtain nu- 
merically the correct perturbation-theory results, and we 
suggest that these tests could be used to benchmark other 
wave-extraction algorithms. 

Finally, we examine the application of the QKT 
scheme to NR simulations of binary-black-hole collisions 
in Sec. VI, considering both the wave zone and the strong 
field regions. We consider a 16 orbit, equal-mass binary- 
black-hole in-spiral and a head-on plunge, merger, and 
ringdown, explicitly illustrating many of the ideas in the 
theoretical discussions of the previous sections. We then 
briefly conclude with a discussion of our results and of 
prospects for the further development of our proposed 
scheme in Sec. VII. 



II. MATHEMATICAL PRELIMINARIES 

In this section, we briefly summarize some important 
properties of Newman-Penrose and orthonormal tetrads 
[Sec. II A], the Weyl curvature tensor [Sec. II B], the 
Lorentz transformations of the Newman-Penrose tetrad 
[Sec. IIC], and the Kinnersley tetrad in Kerr spacetime 
[Sec. II D] . Note that here and throughout this paper, let- 
ters from the front part of the Latin alphabet are used for 
four dimensional coordinate bases, those from the middle 
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part of the Latin alphabet denotes quantities in three di- 
mensional coordinate bases, while Greek indices are used 
for tetrad bases. Bold-face fonts denote vectors and ten- 
sors. 



A. Newman-Penrose and orthonormal tetrads 

Two types of tetrad basis are particularly useful for 
the exploration of generic spacetimes, such as the space- 
times of numerical-relativity simulations of compact- 
binary mergers: i) the Newman-Penrose (NP) tetrad 
basis {e"} = {Z Q , n a , m a , m a }, and ii) an orthonormal 
tetrad {££} = {T a , E 2 ,E 3 ,N a }, which is closely related 
to the NP tetrad as follows. The quantities E 2 and 
E 3 are generally associated with angular variables on 
a closed 2-surface and are related to the complex null 
vector m by E 2 a = y/m(m a ), E% = y/2%(m a ). Here 
3?(m) and 9f(m) denote the real and imaginary parts 
of m, respectively. The timelike vector T and spacelike 
vector N are related to the null vectors I and n by the 
transformations 



r = -^-={T a + N a ) 1 

v 2 



1 



N a ) 



(2.1) 



The metric expressed on the orthonormal basis is the 
Minkowski metric, j al3 = diag{— 1, 1, 1, 1}, while on the 
NP tetrad basis the metric is 
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On the coordinate basis, the components of the metric 
are given by 



ab „ctj9 a b 
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(2.3) 



B. Representations of Weyl curvature tensor 

One aim of this paper is to uniquely fix the NP tetrad 
basis to obtain a set of NP scalars from which an un- 
ambiguous measure of the Weyl curvature (equal to the 
Ricmann curvature in vacuum) can be read off. 

On the NP tetrad, the curvature content of the Weyl 
tensor can be expressed in terms of five complex scalar 
functions 



*0 = 


~C abcd l a m b l c m d 


(2.4) 


*1 = 


-C abcd l a n b l c m d 


(2.5) 


*2 = 


-C abcd l a m b rffn d 


(2.6) 


*3 = 


~C abcd l a n b m c n d 


(2.7) 


*4 = 


-C abcd n a m b n c fn d . 


(2.8) 



An equivalent description of the Weyl curvature can be 
found on the orthonormal frame with associated time- 
like vector T. This is done by defining gravitoelectric 



£ and gravitomagnctic 3 tensors by, respectively, twice 
contracting T with the Weyl tensor and with its Hodge 
dual: 



£ ij — hi" hjC ' a bcdT b T d , 



(2.9) 

By = -^h l a h ] c e abef C ef cd T b T d . (2.10) 

Here h denotes the projection operator onto the local 
spatial hyper-surface orthogonal to T. The normaliza- 
tion for the Levi-Civita tensors is such that 60123 = 1 
and ei23 = 1 hi right-handed orthonormal tetrads and 
spatial triads respectively [sec [16] for a discussion of dif- 
ferent conventions in literature]. These two tensors can 
be combined to obtain a complex tensor 



Si 



iBi. 



(2.11) 



The curvature information contained in Q is exactly the 
same as that contained in the five NP scalars. Recasting 
this information in terms of Q allows us to make use 
of the fact that the £ and 3 tensors describe the tidal 
acceleration and differential frame-dragging to visualize 
the curvature of spacetime [see, e.g., Refs. [16-20]]. 

To make the equivalence between \l/o, ^i, \&2, ^3, ^4 
and Q explicit, we note that the components of the com- 
plex gravitoelectromagnetic tensor expressed on the spa- 
tial triad {E 2 ,E 3 , N} are 



Q 



*2 - 



*0 + *4 



;*o-*4 



*2 + 



+ ^4 



-2* 2 



*l-*3 -*(*l+* 3 ) 

(2.12) 

Furthermore Q is symmetric and trace free (Q\ = 0). 
These results follow from direct substitution of the def- 
inition of the orthogonal basis vectors in terms of the 
NP basis vectors [Eq. (2.1)] and the definition of the 
NP scalars [Eqs. (2.4)-(2.8)] into the definition of Q 
[Eqs. (2.9), (2.10) and (2.11)]. 

Finally, note that for any spacetime in general relativ- 
ity, there are a set of 16 scalar functions or Carminati- 
McLenaghan curvature invariants [21] that can be con- 
structed from polynomial contractions of the Riemann 
tensor. In vacuum, four of these scalars are non- vanishing 
and comprise a complete set of invariants. These four 
scalars can be combined into two complex functions / 
and J and are independent of tetrad choice. In terms 
of the quantities already computed in this section, these 
curvature invariants can be expressed as 
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I=-{£\£\-B\B\)+i{£\B\) 
-^{B\B\B l k -3£\£\B l k ) 

*4 *3 *2 
* 3 * 2 *1 
*2 *1 *0 



(2.13) 
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The invariants / and J play a key role in constructing 
our geometrically motivated coordinate system [Sec. III]. 



C. Lorentz transformations 

There are six transformations of the NP basis vec- 
tors e" that retain the form of the metric given in 
Eq. (2.2). These are the six Lorentz transformations, 
which parametrize the six degrees of tetrad freedom [15]. 
The Lorentz transformations can be decomposed into 
three types depending on which null vector a particular 
transformation leaves unchanged: 



Type I: (I unchanged) 

I — > I , n — > n + am + am 
m — > m + al , m — >• m + al 

Type II: (n unchanged) 



aal 



(2.14) 
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Type III: 

A -2 e 2ie 



A 



-2iQ 



*1 
*3 



* 2 



*2 

(2.19) 



For any algebraically general spacetime, two special 
frame choices exist: the principle null frame (PNF) and 
the transverse frame (TF). The PNF is characterized by 
the property that *4 = = *o; starting from a generic 
tetrad a PNF can be constructed by appropriate Type I 
and Type II Lorentz transformations. The TF is charac- 
terized by the property that *3 = = *i; starting from 
a PNF, a TF can be constructed by additional Type I 
and Type II Lorentz transformations. 

There are three TFs, but only one contains the Kin- 
nersley tetrad in the Kerr limit [10]. In keeping with 
earlier literature [9, 10], we will call this frame the quasi- 
Kinnersley frame (QKF) and the particular tetrad we 
pick out of this frame the quasi-Kinnersley tetrad (QKT). 



Z — >• i + bm + bm + bbn. n —> r, 
m — > m + bn, m — > m + bn 

Type III: (both Z and n unchanged) 



(2.15) 



Z -> A~H, 

iO 

m -> e m, 



n — > An 
m 



e" ie m 



(2.16) 



Here the scalars a and b are complex, while A and G are 
real and can be combined into a single complex number 
A = A^ 1 exp(i0). The rescaling of I and n in Eq. (2.16) 
is called boost freedom, and the phase change of m is 
called the spin freedom. We will follow the convention 
of Ref. [10] and call a set of tetrads related by Type III 
transformations a frame. 

Under the different Lorentz transformations, the NP 
scalars transform as follows: 



• Type I: 

*o 
*i 
*2 

*4 

• Type II: 



*o 
*i 

*2 
*3 
*4 



a*o 
2a* 1 
3a* 2 
4a* 3 



-a 2 * 
-3a 2 *i 
6a 2 * 2 



(2.17) 



a 3 *o 
4a 3 *! 



a 4 * 



D. The Kerr metric and the Kinnersley tetrad 

The no-hair theorems [22, 23] lead us to expect 
all binary-black-hole collisions to ring down to the 
Kerr spacetime after enough time has elapsed. The 
limiting Kerr metric in Boyer-Lindquist coordinates 
(t, r, 9, <P)bl can be expressed as: 



ds = - 1 - 



2Mr 



sm 



, 2 4Marsin 2 6> , , , 
dt 2 dtd(/> 



- a A sin 



A 



dr 2 



(2.20) 



where M and a are the mass and spin of the black hole, 
respectively, and the functions entering the metric are 
defined by 



£ = pp, 



r — la cos ( 



2Mr 



a . 

(2.21) 



For the Kerr spacetime, one tetrad introduced by Kin- 
nersley is particularly conducive for calculation. Among 
other things, on this tetrad the perturbation equations 
in the NP formalism decouple [15, 24]; this feature al- 
lows the perturbation problem to be reduced to the study 
of a single complex scalar (<5*4) that governs the radia- 
tion content of the perturbed spacetime. The Kinnersley 
tetrad expressed on a Boyer-Lindquist coordinate basis 
is given by 



*o 
*i 
* 2 

*3 
*4 



*0 
*1 

*2 + 26*3 
*3 + 6*4 
*4 



46*i + 66 2 * 2 + 46 3 * 3 + 6 4 * 4 



36*2 + 36 2 * 3 + 6 3 * 4 



6 2 * 4 



(2.18) 



l a = I [r 2 +a 2 ,A,0,a] (2.22) 
n a = -L [r 2 + a 2 ,-A,0,a] (2.23) 
[i asin(0),O,l,i csc(0)] (2.24) 
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On the Kinncrsley tetrad, the only non- vanishing NP cur- 
vature scalar is 

* 2 = (2.25) 

In the next section, we explore the behavior of the 
tetrad and curvature quantities defined in this section 
in cases where the physical metric is well understood. 
So doing, we build up some physical intuition that mo- 
tivates our QKT choice, which we then apply to more 
complicated spacetimes, such as those found in numeri- 
cal simulations. 



III. PHYSICAL CONSIDERATIONS FOR 
CHOOSING A TETRAD 

In this section, we introduce several ideas that mo- 
tivate the choice of tetrad and gauge; we will use 
these ideas to explore spacetimes produced by numerical- 
relativity simulations. 

For our purposes, we wish to adopt a tetrad and gauge 
with the following properties (not in order of impor- 
tance): 

1. The tetrad (gauge) reduces to the Kinnersley tetrad 
(Boycr-Lindquist coordinates) when the spacetime 
is a weakly perturbed black hole. 

2. The choice of tetrad and gauge should be indepen- 
dent of the coordinate system, including the slicing 
specified by the time coordinate, used in the NR 
simulation. 

3. To facilitate their real-time computation during a 
NR simulation, all calculation should be local as far 
as possible. 

4. The prescribed use for all computed quantities 
should be valid in strong field regions as well as 
in asymptotic regions of the spacetime. 

5. The choice of tetrad directions should as far as 
possible be tailored to the physical content of the 
spacetime. For example, in asymptotic regions, one 
important direction is that of wave propagation; we 
seek a tetrad that asymptotically is oriented along 
this direction. 

6. To facilitate gravitational- wave extrapolation (from 
the location on the NR simulation's computational 
domain where the waves are extracted to future 
null infinity ^ + ), the falloff with radius of what we 
identify as the radiation field should match that of 
an isolated, radiating system; i.e., it should satisfy 
the expected "peeling properties" . 

We now consider in detail how we may achieve these cri- 
teria in the course of constructing our QKT. 
This section roughly breaks into three parts: 



1. We start [in Sec. Ill A] by motivating the use of 
QKF with a new insight regarding the relation- 
ship between its 1 basis vector and the supcr- 
Poynting vector, which allows it to satisfy crite- 
rion 5. We then review the construction of the QKF 
in Sec. IIIB. 

2. Next, we concentrate on fixing the spin-boost free- 
dom to select the QKT out of the QKF. First of 
all, in Sec. Ill C we discuss several methods for fix- 
ing this freedom that have appeared in literature. 
Then we present our proposal to achieve a global 
and gauge independent fixing [in Sec. HIE] using a 
pair of geometrically motivated coordinates defined 
in Sec. HID. We conclude this part with a brief dis- 
cussion of issues related to the proposed scheme in 
Sees. IIIF and IIIG. 

3. Finally, we discuss [in Sec. IIIH] the conformity of 
the final QKT to criterion 6 and further motivate 
its use. 



A. The TF and wave-propagation direction 

The Kinnersley tetrad [Eqs. (2.22)-(2.24)] is both a 
PNF and a TF [25] [cf. Sec. II C]; this implies that the 
Kerr spacetime is Petrov Type D. Generic non-Type-D 
spacetimes do not have this property: for them no tetrad 
that is both a PNF and a TF exists, so one must decide 
which if either of these properties to preserve. Here, we 
do not want ^4, which plays an important role in the 
perturbation problem, to vanish; therefore, we choose a 
tetrad that is a TF [9-13]. In fact, one particular ad- 
vantage of selecting the TF is its ability to identify the 
direction of wave propagation in the asymptotic region 
[cf. criterion 5]. 

In electromagnetism, a local wave vector that points 
in the normal direction to the surfaces of constant phase 
(wavefronts) can be defined. If the medium through 
which the wave is travelling is isotropic, this direction 
corresponds to the direction of the waves' energy flow, or 
the "wave-propagation direction" , which is determined 
by the direction of the Poynting vector, 

Vi = e ijk E j B k (3.1) 

where the vectors and B k are the electric and mag- 
netic field vectors. In this subsection, we summarize 
the relationship between the QKT and the gravitational 
waves' counterpart to Poynting vector. 

One approach for constructing a geometrically moti- 
vated tetrad follows a suggestion by Szekeres [26] , which 
is to create a gravitational compass out of a number of 
springs. Such a device is sensitive to the spacetime cur- 
vature and can be oriented so that the longitudinal grav- 
itational wave components vanish; mathematically, this 
amounts to reorienting the observer's tetrad so that it 
is a TF, which can be done using Type I and Type II 
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transformations to set v&i = = ^3. We note that 
Chandrasekhar [15] employed the use of a TF for his pro- 
gram of metric reconstruction from a small perturbation 
in curvature 5^ 4 on a background Kerr metric. 

Choosing a TF turns out to orient the tetrad along the 
direction of energy flow, i.e., along the super-Poynting 
vector [27, 28] 

V i = e ijk £ j l B kl , (3.2) 

which defines a spatial direction associated with 
the wave-propagation direction [29]. The super- 
Poynting vector's components in the orthonormal triad 
{Ef,Ef,Ni}, using the explicit form of gravitoelectro- 
magnetic tensor in Eq. (2.12), are 

V E 2 = -P (0, 1) - 3P (1, 2) - 3P (2, 3) - P (3, 4) 
V E 3 = P 1 {0, 1) + 3Pi(l,2) + 3Pi(2, 3) + P(3,4) 

Vn = \ (-|*o| 2 - 2|*i| 2 + 2|* 3 | 2 + |* 4 | 2 ) , (3.3) 

where the functions Po and Pi are defined to be 

P (p,q) = (3.4) 
P 1 (p,q) = R(* p )9f(* 9 )-»(* g )9(* p ). (3.5) 

By transforming to a TF, where ^1=0 = \&3, Eq. (3.3) 
simplifies significantly, becoming 

P = i(|* 4 | 2 -|*o| 2 )iV, (3.6) 

where its direction corresponds to spatial normal direc- 
tion N fixed by our choice of TF and Eq. (2.1), which 
relates TV to the NP tetrad vectors I and n. By selecting 
the TF, we have oriented the tetrad according to the flow 
of energy within the spacetime, achieving criterion 5. We 
believe this is one of the strongest motivating factors for 
making the TF choice. 

B. Computing the quasi-Kinnersley frame on a 
given spacelike hyper-surface 

In this subsection, we review the procedure for con- 
structing the TF that contains the Kinnersley tetrad in 
the Kerr limit. This, as stated before, is named the quasi- 
Kinnersley frame, or QKF. We follow mostly the deriva- 
tion of Rcf. [9]. 

1. A spatial eigenvector problem for the QKF 

Numerical relativity simulations typically split the 4- 
dimensional spacetime to be computed into a set of 3- 
dimensional spatial slices. In the usual 3+1 split, the 
spacetime metric g ab is split into a spatial metric hij, 
lapse a, and shift f3 l according to 

g ab dx a dx b = -a 2 dt 2 + h tJ (dx l + /3 i dt)(dx j + ft dt), (3.7) 



while the Einstein equations in vacuum split into evolu- 
tion equations (for advancing from one slice to the next) 



Rij - ^9ijR = (3.8) 
and constraint equations (satisfied on all slices) 

Rtt - \g TT R = 0, (3.9) 
Rtj - \g Tj R = 0, (3.10) 

where R a b and R are the Ricci tensor and Ricci scalar 
of the spacetime, respectively, the component T is in the 
direction normal to the spatial slice, and the components 
i and j lie within the spatial slice. 

As mentioned in Sec. II, for a given spatial slice with 
future directed unit normal T, the curvature can be ex- 
pressed in terms of the gravitoelectric tensor £ and the 
gravitomagnetic tensor B defined in Eqs. (2.9) and (2.10). 
In terms of the 3+1 quantities typically computed in NR 
codes, provided that the Einstein constraint equations 
are satisfied, the gravitoelectromagnctic tensors in vac- 
uum can be expressed as 

Sa = 3 Rij + KKij - K ik Kj 

Bn = <:'l)i. hi, (3.11) 

where K is the trace of the extrinsic curvature if y , while 
3 Rij and Dfc are the Ricci curvature and connection, re- 
spectively, associated with the spatial metric hij . 

Given the gravitoelectric and gravitomagnetic tensors, 
a powerful tool [16, 17] for visualizing the curvature of 
spacetime is a plot of the "vortex" and "tendex" lines, 
which are the flow lines of the eigenvectors of the grav- 
itoelectromagnetic tensors £ij and Bij. The QKF is 
also related to an eigenvalue problem involving and 
B^, albeit a complex one involving the complex tensor 
Q = £ + iB. Specifically, it was shown in Ref. [9] that 
the QKF can be constructed from the eigenvector a 1 that 
satisfies the eigenvector equation 

Q)ai = -2* 2 a l (3.12) 

where the eigenvalue —2^2 has the value of —2^2 com- 
puted on the QKF. Here and throughout the rest of this 
paper, we adopt the convention of denoting quantities 
associated with a QKF (such as the NP tetrad vector 
I) with an overscript tilde and quantities associated with 
the final tetrad, whose spin-boost degrees of freedom have 
been uniquely fixed (yielding a preferred QKT), with an 
overscript hat (e.g. I). As we will show in greater de- 
tail later in the section, the QKF's Coulomb potential 
^2 can be constructed out of the curvature invariants I 
and J of the spacetime and is invariant under spin-boost 
transformations; therefore, we denote ^2 with a hat to 
indicate it has been fixed to its final value. 
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2. Selecting the correct eigenvalue 



For any symmetric matrix M, the eigenvalues associ- 
ated with the eigenvector problem M^t? = A£ J obey the 
characteristic equation p(X) = where p(X) = det(.M — 
XT) and T is the identity matrix. For a 3 x 3 matrix, the 
characteristic equation becomes 

p(A) = - A 3 + Xhx{M) 

+ \x (tr(M 2 ) - tr 2 (M)) + det(M), (3.13) 



the matrix Q; for a complete derivation, see Ref. [9]. 
The eigenvector corresponding to the eigenvalue —2^>2 
can be expressed as 



if 



(3.17) 



where the real vectors x J and y J are orthogonal with 
respect to the spatial metric hij and their normalization 
obeys the condition 



|5|| 2 -||y|| 2 = i. 



(3.18) 



If Mj = Qj, direct calculation using Eqs. (2.12) and 
(2.13) can verify that tr(Q) = 0, det(Q) = 2J and 
tr(Q 2 ) = 21, which reduces the characteristic polyno- 
mial to 



PeW = -x 3 + xi + 2,j. 



(3.14) 



Here and throughout this section, we will use ||v|| and v ■ 
w to represent norm and inner product of spatial vectors 
under . The vectors x 3 and y 3 can in turn be used to 
define the vectors 



A 1 = 



e ijk x 3 y k 



The solution to this cubic equation can be expressed us- 
ing the speciality index [30] S = 27J 2 /I 3 as 



A = 



I Vs 



(3.15) 



where W(S) = VS - VS~T. There are three solutions 
corresponding to the three transverse frames, but only 
one (namely the QKF) contains the Kinnersley tetrad in 
the Kerr limit [10] (and thus satisfies criterion 1). 

We must now select the correct eigenvalue to define the 
QKF. Only one of the three eigenvalues has an analytic 
expansion around 5=1 (which holds for all Type-D 
spacctimcs, including Kerr [30]). We select this eigen- 
value (which we denote A ) to define the QKF, and so 
—24' 2 = A . For reference, the series expansion of A and 
also the other two eigenvalues A ± around S = 1 is 



A° 



-2*5 



2J 

T 



-3 + §(S- 



1) + 



(3.16) 



X ± ~ ~ 



2. J 
1 



V3 



2 ±l -^ l ~^ S 



In practice, this selection criterion is equivalent to choos- 
ing the eigenvalue with the largest magnitude [10]. 



3. Constructing the QKF tetrad vectors 

We now summarize the necessary results that allow 
the reconstruction of the QKF from the eigenvector of 



1 The fraction on the right of Eq. (3.15) has a three-sheeted Rie- 
mann surface with branch points of order two at S = and 
S = 1, as well as a branch point of order three at 5 = oo. 
The three different eigenvalues arise from the values on the three 
sheets respectively [9], 



/' 



A' + v l + U" 1 X,\ji 



(3.19) 



1 + A • v 

where the normalization condition on <r [Eq. 3.18] ensures 

||A|| = ||^|| = ||5R(/i)|| 2 -||Q(/i)|| 2 = l. (3.20) 

The resulting vectors A, i> and p, turn out to be propor- 
tional to the spatial projections of QKF basis vectors I, 
n and rh respectively. To see this, let the spatial vectors 
be expressed in terms of a spatial triad Ef which is part 
of an orthonormal tetrad with Eq = T a ; then, the 
full QKF tetrad can be constructed as follows: 



s/l - X ■ v 

1-41 

\J\-\-v 
e i0 y] - 



(T a + X i Ef), 
(T a + PE?), 



(3.21) 



+ A ■€> 



V2 y/l-\.v 



Note that the residual spin-boost freedom [cf. Eq. (2.16)] 
has been made explicit in Eq. (3.21) by means of the 
parameters A and 8 (which have yet to be determined) . 

Also note that the equation for rh above must be mod- 
ified if the normal to the spatial slice T lies in the plane 
spanned by I and n, since in this special case the vectors 
A and v turn out not to be independent of each other (as 
is true generally) but are instead related by A 4 = —v % . 
It is unclear whether such a slicing can be found for any 
spacetime, but once found, it is closely associated with a 
TF. In this case the vector p, is undefined and rh should 
be constructed from any real unit vector f in the spatial 
2-plane orthogonal to A and T according to 



iO 



—(r i + i j*\ j r k )E?. 



(3.22) 
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Because the spatial eigenvector problem (3.12) can be 
solved point-wise, the construction of the QKF is a lo- 
cal procedure and satisfies criterion 3. Furthermore, the 
procedure can be applied in the strong field region [cf. 
criterion 4], although the physical interpretation is only 
clear if the tetrad can be smoothly extended from there 
to infinity. By choosing our tetrad to be a QKF, we have 
used up four of the six possible degrees of tetrad freedom 
and have uniquely fixed the directions associated with the 
real null vectors I and h. We will address the remaining 
spin-boost freedom in the next three subsections. 

C. The spin-boost tetrad freedom 

After electing to work in the QKF, the residual tetrad 
freedom is restricted to a Type III spin-boost transforma- 
tion [Eqs. (2.16) and (3.21)]. As seen in Eq. (2.19), the 
boost transformation affects the magnitude of Vt^, while 
the spin transformation modifies the phase of 

To gain some insight into what the spin-boost transfor- 
mations do physically, consider a congruence of observers 
whose world lines are the integral curves of the T field in 
Eq. (2.1). For these observers, a spin transformation of 
the tetrad mixes up the two polarizations of gravitational 
wave by the induced phase rotation 2 ; in practice, this 
rotation occurs because the observers are rotating the 
orientation of their coordinates and thus redefining what 
they consider to be the latitudinal and longitudinal direc- 
tions. Similarly, the boost transformation in Eq. (2.16) 
alters the velocity with which these observers move along 
the direction of wave-propagation, causing the gravita- 
tional wave they observe to be redshiftcd or blueshifted. 

In order to identify the gravitational wave and curva- 
ture content contained in "J 4 in an unambiguous man- 
ner, we need to provide a prescription for fixing A and 
throughout the spacetime. Note that A and v constructed 
in Eq. (3.19) are dependent on the choice of slicing; thus 
simply setting A and in Eq. (3.21) to particular values 
does not select a tetrad in a slicing independent manner. 
Fixing these parameters but altering the slicing will lead 
to different tetrads in the same frame (the QKF), thus 
when we leave A and undetermined, the frame as a 
whole that we obtain from Eq. (3.21) is slicing indepen- 
dent. 

One example of fixing the spin-boost freedom in a 
gauge independent way often used in mathematical anal- 
ysis is selecting the so-called canonical transverse tetrad 
(CTT) [25], which is defined by the condition that 

*i = = * 3 and * = * 4 . (3.23) 

The CTT has the property that the super-Poynting vec- 
tor given in Eq (3.6) has vanishing magnitude; in this 



tetrad, the observers are co-moving with local wave- 
front in the asymptotic region and consequently measure 
HP 1 1 = 0. Since no physical observer can travel at the 
speed of light and co-move with the wavefront, we re- 
quire a more physically motivated prescription for fixing 
the spin-boost freedom. 

Several approaches for providing such a physically mo- 
tivated prescription have been suggested. A common ap- 
proach is to impose conditions on spin coefficients (such 
as e = [8]). The Kinnersley tetrad for the Kerr metric 
has spin coefficients 3 that obey n = <7 = \ = v = e = 0. 
The meaning of some of these coefficients can be gleaned 
from the equations governing how the tetrad evolves 
along the I direction, namely [15] 

1% = 2K(e)Z a -rem a -rem Q , (3.24) 
l b m? b = 2i%(e)m a +Wl a - nn a . (3.25) 

If re = for example, the null vector I is tangent to a 
geodesic and further if 9?(e) = this geodesic is affincly 
parameterized. 

Note that choosing I to be geodesic or re = is not 
necessarily consistent with choosing to work in a TF, al- 
though these conditions are consistent in the Kerr limit. 
In a TF, the only freedom available to set the spin co- 
efficients to zero is the spin-boost transformation. Since 
re transforms as re — > A~ 2 e lB K under Eq. (2.16), the spin 
coefficient re cannot be set to zero. The spin coefficient 
e, on the other hand, transforms as 

e -> A-h - ^A- 2 l a V a A + ^A-H a V a e, (3.26) 

and can be made to vanish by suitably chosen A and 0. 
Equations (3.24) and (3.25), indicate that the condition 
e = can be used to fix the scaling of I as well as the 
phase of m. Setting e — can therefore be used as a 
means of fixing the spin-boost freedom, but this choice 
has the disadvantage that Eq. (3.26) must be solved in 
order to obtain A and 0, which can be expensive numer- 
ically. 

In the following subsections, we present an alternative 
method of fixing the spin-boost freedom by constructing 
a coordinate system based on the curvature invariants. 
Differentials of these new coordinates are then used to 
set the scale or fix the spin degree of freedom of the final 
QKT. This method avoids the need to solve differential 
equations by directly imposing local conditions of the the 
tetrad basis vectors. 



D. A geometrically motivated coordinate system 

In this paper, we fix the spin-boost freedom by exploit- 
ing the curvature invariant "J 2 [identified in Eqs. (3.12) 



2 Recall that for plane waves on Minkowski background, we have 
\I/4 = — h+ + ihx, where h is the metric perturbation. 



3 For how the spin coefficients [which are complex scalars] are de- 
fined in terms of the null tetrad, see e.g. Eq. (1.286) of Rcf. [15]. 



FIG. 1: Properties of the (f , 9) coordinates constructed from the Coulomb potential in the QKF. (a) The equatorial plane of a 
Kerr spacetime in a Kerr-Schild slicing with contours of constant Boyer-Lindquist radius f at equal increments. The inset zooms 
in around the event horizon (indicated by a transparent black disk). The f contour increments in the inset, while still uniform, 
are smaller than in the main figure, and the thick contour line coinciding with the event horizon matches the Boyer-Lindquist 
radius f+ in Eq. (3.30). (b) Surfaces of constant latitudinal coordinate 9 for the Kerr-Schild slicing. 



and (3.16) and computed using Eq. (3.15)] to define geo- 
metrically motivated and unambiguous radial and latitu- 
dinal coordinates. The quantity ^2 can be interpreted as 
the Coulomb potential experienced by an observer [26], 
and all observers in a QKF agree on its value. Our pre- 
scription for fixing the spin-boost freedom is to effectively 
tether our observers to a fixed position with respect to 
the coordinates associated with the instantaneous back- 
ground Coulomb potential they experience. By doing 
this, we choose "stationary" observers that watch grav- 
itational waves pass, in contrast to the CTT observers 
(Sec. IIIC) that co-move with the waves. In Kerr limit, 
our choice amounts to selecting a set of stationary ob- 
servers associated with the Boyer-Lindquist coordinate 
system. 

To illustrate this idea more fully, note that when we 
work within the QKF, the complex gravitoelectromag- 
netic tensor from Eq. (2.12) reduces to 



Q 



*2 - (#0 + *4)/2 i(#o-#4)/2 
i(*o-*4)/2 * 2 + (*o + * 4 )/2 

-2*2. 

(3.27) 

making TV an eigenvector. Of particular interest is the 
component 



NN 



^ NN ' 



(3.28) 



As illustrated in detail in [17] and particularly in Sec. IV 
A of Ref. [16], within the context of vortexes and ten- 
dexes, £^ N measures tidal acceleration and B^^ the dif- 
ferential frame-dragging experienced by a person whose 
body is aligned along the radial TV eigenvector. The 
frame dragging induced by the angular momentum of 
the source implies a latitudinal coordinate, and the ra- 
dial tidal acceleration implies a radial coordinate. The 
Coulomb potential ^2 thus contains information about a 



pair of geometrically motivated coordinates f and 9. 

To relate the Coulomb potential ^2 to the geometric 
coordinates f and 9 in a meaningful way that reduces to 
the Boyer-Lindquist coordinates in the Kerr limit (thus 
satisfying criterion 1), we make use of expressions for 
the Kerr spacetime [Eqs. (2.25) and (2.21)] to define the 
coordinates. In other words, we define f and 9 using the 
complex equation 



p = f — ia cos(0) 



M 



1/3 



(3.29) 



where M and a are real constants that become just the 
mass and spin of the central black hole in the Kerr limit. 
A discussion regarding these parameters in dynamical 
simulations follows in Sec. IIIF. Recall that the Coulomb 
potential ^2 can be constructed directly from the curva- 
ture invariants / and J of the spacetime; the construction 
of the coordinates out of curvature invariants makes them 
slicing or gauge independent, thus satisfy criterion 2. 

Figures 1 and 2 explore some properties of f and 9. 
The first property is the ability to recover the Boyer- 
Lindquist radial and latitudinal coordinates from a Kerr 
spacetime expressed in any slicing. A particular exam- 
ple using Kerr-Schild slicing is shown in Fig. 1, where we 
plot the contours of f and 9 under Kerr-Schild spatial 
coordinates (r, 9, (f>). The resulting figures show that the 
coordinate transformations between (f,6) and (r,6) (un- 
like those for Boyer-Lindquist t and </>) do not become 
singular at the event horizon [cf. criterion 4], which co- 
incide with the contour of 



f + =M 



M 2 - a 2 



(3.30) 



as expected. The (f,9) coordinate system for a dynami- 
cal simulation of two equal-mass, nonspinning black holes 
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FIG. 2: (a) Snapshot surfaces of constant r for an equal-mass nonspinning binary merger simulation taken during the inspiraling 
phase. Far away from the black holes, the contours represent those expected from a monopole moment. When moving closer 
to the black holes, higher order multipoles become important, (b) Constant 6 surfaces for the same simulation as in (a), shows 
a "spiral-staircase" pattern generated by rotating deformed cones as discussed in greater detail in Sec. VIA. 



during their inspiral phase is shown in Fig. 2. The peanut 
shaped features in panel (a) makes apparent the fact that 
the coordinate system is adjusting to the intrinsic geom- 
etry of the simulation. The cones of constant angular 
coordinate 6 display a wavy feature when compared to 
the simulation coordinate 9. This feature and its origin 
will be discussed in greater detail in Sec. VIA, where we 
explore the binary simulation in more detail. 



E. Fixing the spin-boost degrees of freedom 

The previous subsection provides us with an unam- 
biguous and geometrically motivated set of radial and lat- 
itudinal coordinates that are valid throughout the space- 
time and that are independent of the choice of slicing. 
Our strategy for fixing the last two degrees of tetrad free- 
dom is to require that the tetrad frames can be associ- 
ated with observers that are in some sense "stationary" 
with respect to our geometrically motivated coordinates 
while also requiring that the selected tetrad reduces to 
the Kinnersley tetrad in the Type-D limit. 

To achieve this construction (and thus to provide a 
global prescription for fixing the spin-boost freedom), 
note that df provides a measuring rod in the radial di- 
rection, relative to the wavefront, against which the scale 
of the radial component of I can be fixed. Similarly d9 
provides a transverse direction which can be used to fix 
the phase of rh. Let us now begin with any tetrad in the 
QKF {l,n,fh,rh}, constructed according to Eq. (3.21). 
The prescription we use to fix the parameters A and 
associated with the spin-boost degrees of freedom to ob- 
tain the final QKT {l,n,rh,rh} is to require that the 
final tetrad obeys 



arg 



(df)J a = h 
(ddj m a = arg[p] 



(3.31) 
(3.32) 



Note that these conditions are exactly the conditions sat- 
isfied by the Kinnersley tetrad in Eq. (2.22) and (2.24) 
except that the Boyer-Lindquist coordinate has been re- 
placed by its corresponding geometrically constructed 
counterpart introduced in Sec. Ill D. The reduction to 
the Kinnersley tetrad in the Type-D limit is thus triv- 
ial [cf. criterion 1]. Furthermore, conditions (3.31) and 
(3.32) contain only local differentiation and algebraic cal- 
culations and thus obey criterion 3. They also inherit 
gauge independence from the QKF and the geometric 
coordinates, thus satisfy criterion 2. 

It turns out that the final QKT can be constructed by 
starting with a tetrad in the QKF with A = 1 and 9 = 
in Eq (3.21), computing the quantities 



A 



(df)J a , 



O = — are 



(^) a m a ] +arg[ / S], 



(3.33) 
(3.34) 



and then substituting these values back into Eq. (3.21) 
to obtain the final tetrad. Our fictitious observers have 
now oriented and scaled their tetrads according to the 
Coulomb potential they experience by observing the lo- 
cal changes in tidal acceleration and differential frame 
dragging. 



F. The effect of a and M on the tetrad choice 

In the definition of the geometric coordinates (f, 6) 
in Eq. (3.29), two constants M and d corresponding to 
the mass and spin of a Kerr black hole in the Type-D 
limit entered our prescription. We now clarify their in- 
fluence on the final computed quantities of $4 and the 
constructed tetrad. 

First, we observe that the spin a docs not affect the 
spin parameter O in expression (3.34) and can be left 
undetermined, since only the direction of d6 is required 
to determine the argument of its inner product with rh. 
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The final computed quantities are however dependent 
on the value of M, which enters as a constant factor scal- 
ing the boost parameter A. The computed ^4 is simply 
rescaled by a constant scaling factor if the value of M is 
changed. This allows one to compute all quantities real 
time during the simulation with (say) M = 1 and to a 
posteriori rescale the results once the final mass of the 
remnant black hole is known. 



G. The remaining gauge freedom 

Using the appropriate combination of the curvature 
invariants [Sec. Ill D] to prescribe radial and latitudinal 
coordinates (f , 9) fixes two of the four degrees of gauge 
freedom, while the choice of a TF [Sec. IIIB] and the 
subsequent fixing of the spin-boost freedom [Sec. Ill E] 
removes all six degrees of tetrad freedom. What remains 
is to fix the final two degrees of gauge freedom: the slicing 
(or time coordinate t) and the azimuthal coordinate 4>. 

For a given slicing, "far enough" from the strong field 
region, surfaces of constant r and 9 intersect in a circle. 
This can be seen graphically in Fig. 2 by superimposing 
plot (a) and (b) and taking "far enough" to mean the re- 
gion where the mass monopole and current dipole are the 
dominant terms in the Coulomb background. The pre- 
scription of the azimuthal coordinate q> is then as simple 
as requiring that given a specific (as yet undetermined) 
starting point, the proper distance increments dip along 
the circle remains constant. 

Fixing the time slicing requires more finesse. One 
method of specifying the time slicing indirectly is by 
means of a congruence of outward propagating affinely 
parameterized null geodesies [see Sec. IIIH2 below for a 
suggested congruence] starting from a fixed radius f; the 
affine parameter r is then used as a coordinate. This ap- 
proach is particularly suited to the task of wave extrac- 
tion where the quantities computed should exhibit the 
scaling laws predicted by the peeling property [31, 32]. 

The prescriptions given above contain residual free- 
dom. Fixing them is beyond the scope of our current 
work. In this paper, wherever needed, we simply use the 
coordinate time in the simulation and the simulation's 
azimuthal coordinate. 



H. The peeling theorem 

1. Peeling in Newman-Penrose scalars 

In this section, we consider the peeling property, which 
describes the way in which, for an isolated gravitating 
system that is asymptotically flat, the components of the 
curvature tensor fall off as one moves farther away from 
the source of the emitted gravitational radiation. At suf- 
ficiently large distances, only Type N radiation is notice- 
able; the limiting Type N radiation can be identified as 



the gravitational-wave (GW) content of the spacetime 
(typically denoted as ^4 on an affinely parameterized 
out- going geodesic null tetrad). [Note that gravitational 
radiation is only rigorously defined at future null infinity 
(denoted J f+ )-] A caricature of this behavior is given in 
Fig. 3. 

Here we review the usual derivation of the peeling 
property [25, 31-34], commenting on some of the proper- 
ties of the QKT within this context; an alternate deriva- 
tion of the the peeling property using spinor notation can 
be found in [34] . The basic idea of the usual derivation 
is to introduce a new 'unphysical' metric ds that is con- 
formally related to the physical metric ds by ds = £1 ds. 
The metric ds is finite and well defined where the physi- 
cal metric blows up (points on ^ + are infinitely distant 
from their neighbors [25]) and allows us to explore the 
properties of the spacetime at J? + or at conformal null 
infinity, where O — >■ 0. All quantities associated with the 
conformal metric ds will be denoted with an acute (e.g. 
ds). 

The relationship between metric tensors can be ex- 
pressed as 

9ab = tfg ab , g ab = rrV", (3.35) 

and the topology at J + is S 2 x R. Now let l a be tangent 
to an affinely parameterized out-going null geodesic on 
the real spacetime, with an affine parameter r such that 
l a V a T = 1. Then let l a be tangent to an affinely param- 
eterized geodesic in the conformally related spacetime 
with affine parameter f. Note that if we take l a = fl 2 l a , 
then the geodesic equation in physical spacetime implies 
its counterpart in the conformal spacetime [25]; further- 
more, if we choose n a = n a at J^+j then we have that 
the direction of {n a = n a )\j?+ does not depend on the 
geodesic and is tangent to [25]. 

Substituting these choices into the expressions for the 
metric [Eq. (2.3)] and subsequently into Eq. (3.35) we 
have that at ^ + the conformal tetrad relates to the phys- 
ical tetrad by means of the expressions 

l a = n 2 l a , m a = nm a , n a = n a . (3.36) 

Departing from J? + by moving into the manifold, differ- 
ences in parallel transport in the physical and conformal 
manifolds lead to higher order terms in the m and n 
equations (see Eqs. (9.7.30) and (9.7.31) in Ref. [25]). 
By comparing the affine parameter on the two manifolds 
along a geodesic and imposing Einstein's vacuum field 
equations, we can show that in general dr = Q,~ 2 di and 
that for large affine parameter r or small conformal affine 
parameter f we have [25] 

f = -A-^t' 1 + D nT~ n , (3.37) 

n=2 

T = -A- 2 T- 1 + J2 C nT n , (3.38) 

= A- x t- x + ^£„t-™, (3.39) 

n=2 
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Q, = -Af- y £ d A n T n , (3.40) 

n— 3 

where A n , C n , D n . E n are constants and A = — 4§|f2_».o 
is a non-zero constant. Any quantity 0... that is C h con- 
tinuous at ^ + can be expressed in terms of a series ex- 
pansion about as follows 

h 

9... = Y / r n i -- ) +o(r h ) 

h 

= ^T- n d { " ) +o(T- h ). (3.41) 
n=0 

Since the Wcyl tensor is conformally invariant, C a bcd = 
<Kd> or 

Cabcd = ^ 2 C a bcd, (3.42) 

all the relevant quantities can be computed on the confor- 
mal manifold where the metric is finite and well behaved, 
and then interpreted on the physical manifold where 
the metric quantities may have diverged. At in an 
asymptotically fiat spacetime, the Weyl tensor C a bcd van- 
ishes and the dynamics of the gravitational field as one 
approaches ,y + can be described using a tensor K a bcd, 
where 

Cabcd = QKabcd (3.43) 

and the components of K expressed on the tetrad basis 
{I, n, rh, m} admit expansions in the form of Eq. (3.41). 

The peeling-off property of the Weyl scalars naturally 
arises when one expresses the quantities related to K 
in terms of the physical metric and the tetrad basis 
{I, n, m, m}. Let us take a detailed look at ^4: anal- 
ogous to the definition of "J 4 in Eq (2.8), let 

*4 = -K a bcdn a m'n c rn i (3.44) 

The fact that K is regular as we approach J> + implies 
that ^4 admits a series expansion of the form 

*4 = $>~"*i n) ' (3.45) 

n=0 

where in particular = ^^\jt+. Similar expansions 
can be found for 4>i,i = 0,1,2,3. At J^ + , the physical 
* 4 [defined by Eq. (2.8)] is related to ^4 by 

* 4 = -(n- 2 c abcd )(n a )(nm\n c )(nm) = fi* 4 , 

(3.46) 

where we have used Eqs. (3.36), (3.42) and (3.43). By a 
similar argument as used for ^4, the differing powers of 
CI appearing in Eq. (3.36) result in a hierarchy being set 
up where 

&i = ft 5 ~%. (3.47) 



This expression is merely a product of the series in 
Eq. (3.39) and (3.41). Rcsumming the product of series 
implies that the physical Weyl scalars along an affincly 
parameterized out-going null geodesic can be expressed 
as 

*i = r'- 5 T ' n 4 n) (3-48) 
n=0 

where tp^ are constant along the geodesic. 



2. Peeling in principal null directions 

Note that the peeling property is not a function of 
which geodesic is chosen (provided that the geodesic 
strikes and is affincly parameterized); on the con- 
trary, it is a feature of the spacetime curvature and the 
distribution of principal null directions (PNDs) as one 
approaches J^ + . This feature is illustrated graphically in 
Fig. 3 (a): as one moves in toward the source from 
along a null geodesic, the PNDs "peel off" away from the 
geodesic direction [34]. 

Let us now quantify this behavior more precisely. 
Starting from the I vector associated with the out-going 
null geodesic, perform a Type II Lorentz transformation, 
so from Eqs. (2.15) and (2.18) we have that the four prin- 
ciple null direction (PNDs) can be expressed as: 

k = I + bm + bm + bbn, (3.49) 

where b takes on the values of the four roots of the com- 
plex equation 

*o + 46*i + 66 2 * 2 + 4fe 3 * 3 + 6 4 * 4 = 0. (3.50) 

From Eq (3.49) it becomes apparent that the magnitude 
of b determines the extent to which the PNDs depart 
from the null vector I since k a l a — —bb. By making the 
identification proposed in [10] between a pair of spherical 
coordinates (9, 4>) and the boost 6, 

& (4) = cot (§) e ^> i€ {1,2,3,4}, (3.51) 

we can graphically demonstrate the motion of the PNDs 
by plotting the four roots on the anti-celestial sphere as 
shown in Fig. 3 (b). (The anti-celestial sphere can be 
thought of as the space of all possible directions asso- 
ciated with out-going null rays.) If 9i — 7r, then the 
magnitude of the boost bu\ vanishes and k = I is a PND; 
on the other hand, if 9{ = then k oc n. 

Asymptotically, where the Weyl scalars admit power 
series expansions such as Eq. (3.48), we can obtain the 
dominant behavior of b by setting 

b = T~ n b^ (3.52) 

71 = 
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FIG. 3: (a) A pictorial representation of the peeling property as bunching of principal null directions (PNDs) [25], with the inset 
showing a top-down view, (b) The relationship between the PNDs and the quasi-Kinnersley tetrad (QKT). Points A,B,C,D 
correspond to PNDs on the anti-celestial sphere, and are arranged as the vertexes of a tetrahedron. The anti-celestial sphere 
can be thought of as a spatial slice of the future null cone, where each point on the sphere represents a null direction. The 
line EF linking the mid-points of a pair of opposite edges strike the anti-celestial sphere at opposite poles which correspond to 
the null direction I associated with the QKT at point (F) and to the null direction n associated with the QKT at point (E) 
[cf. Fig. 8-5 in Ref. [25]]. (c) The four principal null directions recorded during a head-on numerical simulation [described in 
Sec. VI B] are represented as points (in four different colors) on the anti-celestial sphere. We begin integrating a null geodesic 
in the I direction and then compute the PNDs at discrete intervals along that geodesic. Darker colored points correspond to 
values farther along the geodesic (farther removed from source region). For cleaner visualization, the angular coordinates on 
the anti-celestial sphere in this figure are simply those of the simulation coordinates and not the abstract ones in Eq. (3.51). We 
nevertheless see that the PNDs are distributed in a pairwise symmetric manner relative to tangent £ of the geodesic (denoted 
by the black radial line). Two of the PNDs stay close to £, whose close-ups are shown in the framed inset. The other two 
demonstrate a clear motion toward I, where arrows indicate progress along the null geodesic. The numerical findings are thus 
consistent with the bunching behavior depicted in panel (a). 



and substituting this expression into Eq. (3.50). We then 
have that b^> = and M 1 ) can be found by finding the 
four roots of the equation 

4°) + 4&«v4°> + e(V l >) 2 v4 0) 

+ 4(fo( 1 >)%f + (V 1 >)%f =0. (3.53) 

Further higher order terms become more complicated and 
involve mixtures of higher order terms in the expansions 
of the Wcyl tensor components. 

The leading order coefficients tpf^ in Eq. (3.48) are in- 
dependent of the choice of geodesic path, while higher 
order terms ^| with n > arc path or geodesic- 
dependent, which implies in turn that the are 
geodesic-dependent. This path dependence suggests the 
possible existence of an optimal null trajectory along 
which the series converges most rapidly and from which 
the GW content can be most effectively extracted. One 
approach to finding the optimal trajectory is to minimize 
the higher order terms, ip^ (n > 0), achieving a rapidly 
converging series. Possibly the most rigorous method of 
ensuring rapid convergence would be to identify the Kin- 
nersley tetrad and thus the wave propagation direction 
at ^ + and then to integrate backward in time, but such 
a strategy cannot be executed real time during a numeri- 
cal simulation. Instead, the method advocated here is to 



align the initial geodesic direction with the wave propa- 
gation direction in the computational domain and then 
to integrate forward in time. This direction can be iden- 
tified in a slicing independent way by I in the QKT as 
was shown in Sec. Ill A. In Sec. VI, we will demonstrate 
numerically the rapid convergence rate that results from 
this approach. 

Choosing the QKT I as the initial direction is further 
justified by considering the manner with which PNDs 
converge onto the outgoing gcodesic's tangent direction. 
In the QKT = = \&3, which greatly reduces the 
complexity of Eq. (3.50). The transformation from I to 
PND takes the simplified form 

b 2 = -J- (^-3*2 ± - *4*o^) • (3.54) 

The four roots now occur in pairs and can be parameter- 
ized using only two angles. 

S (1) =cotf|)e* 6 (2) = cot (| ) e*+* 

S (3) =cot(|)e* S w =cot^e*+«" 

(3.55) 

The out-going null direction I of the QKT thus finds it- 
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self in the center of the four PNDs due to the added sym- 
metry imposed by the QKT. This situation is depicted 
graphically in Fig. 3 (b). By initially selecting a QKT di- 
rection in the interior of the computational domain from 
which to shoot the geodesies to infinity, we impose an 
additional symmetry on the manner in which the PNDs 
approach the gcodesic's tangent initially, hoping that this 
additional symmetry is maintained as the geodesic ap- 
proaches to ensure the clean pair wise convergence of 
the PNDs to the geodesic's tangent. 

Once the geodesic is shot off in the I direction, there 
is nothing to ensure that it remains in the QK out-going 
null direction. In practice, however, the QK property 
appears to be maintained to a high degree of accuracy, 
as is indicated by the symmetric pairwise convergence 
of the PNDs onto the null geodesic shown in Fig. 3 (c). 
For this plot the angle between the QKT direction of 
I and the tangent t to the geodesic remains less than 
4.2 x lCTSr. 



3. Peeling of QKT quantities 

We close this section on the peeling property by re- 
visiting the geometrically motivated coordinate system 
(introduced in Sec. Ill D) in the asymptotic region. The 
curvature invariants / and J (and thus $2) can be con- 
structed using the series expressions Eq. (3.48). The 
dominant behavior of the curvature invariants arc 



~ 3 ^ (3.56) 



[see Eq. (2.13)] where the quantities with a super- 
script are constant along the geodesic. Assigning the 
radial coordinate using Eq (3.29) sets 



t3£ 



1/3 



(3.57) 



The peeling property states that the PNDs converge onto 
the out-going geodesic direction £. Since each pair of 
PNDs arc equidistant from the QKT I, this implies that 
I approaches the £ direction. The asymptotic relationship 
between f and r given in Eq. (3.57), together with the 
condition Eq. (3.31) that we use to fix the boost freedom 
of the QKF, implies that I not only asymptotes to the di- 
rection of £, it is also affinely parameterized in this limit. 
The geometrically constructed f asymptotically denotes 
the spherical wavefronts of light-rays approaching . 

Lastly, we underscore the fact that using the QKT has 
the advantage of identifying a unique affine parameteri- 
zation of the geodesic as it approaches . The prescrip- 
tion given in Eq. (3.31) for fixing the boost freedom of 
the QKT has used the geometry of the spacetime implicit 
in the Coulomb potential to fix the parameterization of 
I in a global manner, removing the freedom to choose 
a different affine parameter through the transformation 
t — > At. These ideas will be revisited in greater de- 
tail when we look at extrapolation in the context of the 
numerical simulations in Sec. IV B. 



IV. NUMERICAL IMPLEMENTATION 

In this section, we detail the numerical implementation 
of the analytic ideas mentioned in the previous sections 
using the Spectral Einstein Code (SpEC). A description 
of SpEC and the methods it uses are given in Ref. [35] 
and the references therein. 



A. Constructing the QKT 

We construct the QKT in a numerical simulation by 
first constructing an orthonormal tetrad adapted to the 
simulation's coordinate choice and then the orthonor- 
mal tetrad's null counterpart {I, n, m, rrf] and the 
associated NP scalars 4^. In order to find a QKF 
{I, n, m, rh}, the construction described in Sec. IIIB 
can be used; alternatively, the appropriate Type I and 
Type II transformations [Eqs. (2.14) and (2.15)] to the 
QKF can be found. Finally, we construct the geomet- 
rically motivated coordinate system (r, 9) described in 
Sec. Ill D, and we use these coordinates to fix the remain- 
ing Type III tetrad freedom to obtain the QKT. 



1. Implementing a coordinate tetrad 

Specifically, we begin our construction by noting that 
the SpEC code stores the spacetime metric g a b on a 
Cartesian coordinate basis {x a } = {t, x, y, z}. (Note 
that henceforth the index refers to the time coordinate.) 
We can also define a set of related spherical coordinates 
{t, r, 9, 4>} by using the standard definitions 



r sin cos < 



r sm a sin ( 



rcos6». (4.1) 



We further define the time-like unit normal to the spatial 
slicing and radially outward-pointing vector as 



ra _ aa 

rpa _ "0 P yya _ 



(4.2) 



respectively, where a is the lapse and (3 a is the shift, and 
r is the spatial location vector. Inserting these orthonor- 
mal vectors into Eq. (2.1) yields 1 and n, two legs of the 
null tetrad tied to the simulation's coordinates. 

We next construct the remaining two tetrad legs 
{m, m}, ensuring that the normalization conditions of 
Eq. (2.2) are satisfied. In other words, we seek to con- 
struct the null vector m = 1/^2(E 2 + iE 3 ) where E 2 
and E 3 are orthogonal to T and N and to each other 
and obey the normalization condition 



\E 2 



\E 



1. 



(4.3) 



Our construction begins by computing the vectors 

19 n 13 
— : > F = , 

r sin 9 d4> ' r 89 ' 



K 



(4.4) 
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where 6, <j> are spherical coordinates defined in Eq. (4.1). 
Then, we ensure orthogonality by means of the Grams- 
Schmidt-likc construction 



(F) a = F a + F b l b n a + F b n b l a , 



(4.5) 



rescaling appropriately to obtain the correct normaliza- 
tion as follows: 



(E 



2\a 



pa 



(4.6) 



Similarly, for the final tetrad leg, we construct the or- 
thogonal vector 

(K) a = K a + K b l b n a + K b n b l a - K b E b (E 2 ) a , (4.7) 

normalizing it as follows: 

K a 



(£ d ) a = 



(4.8) 



2. Obtaining a tetrad in the QKF 

Given the orthonormal coordinate tetrad 
{T a , N a , (E 2 ) a , (E 3 ) a } : we next construct a tetrad 
{I, n, m, m} in the QKF by using the results of 
Sec. IIIB, in particular Eqs. (3.12), (3.17), (3.19) and 
(3.21). We can alternatively construct a QKF tetrad 
by explicitly rotating our initial coordinate tetrad into 
a transverse one via Type I and II transformations 
[Eqs. (2.14) and (2.15)]. We have implemented both 
constructions numerically and verified that they agree; 
in the remainder of this subsubsection, we discuss details 
of each implementation in turn. 

The hyper-surface approach of Sec. Ill B requires us to 
solve the complex eigenvector problem in Eq. (3.12), with 
Q calculated cither from Eq. (3.11) or from Eq. (2.12). 
Using Eq. (2.12), the eigenvector problem can be solved 
analytically. After computing the desired eigenvalue 
A = — 2^2, which is the root of Eq (3.15) that admits 
the expansion (3.16) (in practice, it suffices to select the 
eigenvalue with the largest norm as suggested by Beetle 
ct al. [9]), the corresponding un- normalized eigenvector 
£ of matrix (2.12) is 



t E 2 =2A 2 - A (¥ + ¥ 4 - 2* 2 ) 



+ 2 (¥ x + ¥ 3 ) - ¥ 2 (¥ + ¥ 4 + 2¥ 2 ) , 

S B 3 =i [A (¥ - *4) + 2 (¥§ - ¥? + ¥ 2 (¥ - ¥4))] , 

Stv =2 [(¥ + ¥ 2 ) ¥3 + A (¥1 - ¥3) - ¥l (¥ 2 + ¥4)] • 

(4.9) 

where the ¥j values arc those extracted on the coordinate 
tetrad. (Note that this formula fails when ¥1 = = ¥3, 
but in this case the coordinate tetrad is already in the 



QKF.) To normalize X into <x that satisfies Eq. (3.18), 
we multiply it with a suitable complex number, namely 



J_ /|o|^ VW+J_ + _\/2j< 

V2 V < 



where 



E" 

/3 



7V/3 + 7 



X a + iY a , a = X a Y a , 
||X|| 2 -||11| 2 , 1=VW 



4a 2 



S (4.10) 



(4.11) 
(4.12) 



Alternatively, we can construct the QKF using the 
Type I and II transformations applied to the coordi- 
nate frame as follows. Starting from a general Petrov 
Type I spacetime with five non-vanishing Weyl scalars, 
we perform a Type I rotation, introducing a parameter 
a, followed by a Type II rotation that introduces a pa- 
rameter b. These parameters can then be chosen to set 
¥1 = ¥3 = by solving the resulting system of two 
equations for the two parameters a and b. Reference [10] 
shows that the appropriate choice of parameters can be 
found by defining the intermediate quantities 

H = * *2 - *?, G = ^l^ 3 - 3*0*1*2 + 2*J 



(4.13) 



and then setting 



_ G± v /G 2 + (*oA-2g) 2 (g + * A) 



* n A - 2H 



¥3 + 3a¥ 2 + 3a 2 ¥i + a 3 ¥ 



(4.14) 



¥4 + 4a¥ 3 + 6a 2 ¥ 2 + 4a ci ¥i + a 4 ¥ 

(4.15) 

Note that this prescription becomes ill defined when ¥0 
on the the initial tetrad approaches zero or when ¥oA — 
2H = 0, making it difficult to find a by solving Eq. (4.14); 
this problem is easily resolved by first applying a Type 
II transformation that takes the initial tetrad into one 
in which these pathologies do not arise. Furthermore, 
we have two possible solutions for a resulting from the 
freedom to interchange the I and n legs associated with 
the transverse frame^ the convention we use is to choose 
the root that gives (l — n) a l a > 0, i.e. we choose I to be 
outgoing in the simulation coordinates. 



3. Obtaining the quasi- Kinnersley tetrad from the 
geometric coordinates 

With a QKF in hand, we next seek to specialize to 
the particular QKT described in Sec. HIE, where we 
use geometrically motivated coordinates (f, 6) given by 
Eq. (3.29) to fix the final Type III degrees of freedom. In 
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order to fix these freedom using Eqs. (3.31) and (3.32), 
we must calculate the one-forms df and d9. We com- 
pute the spatial derivatives spectrally, and we compute 
the time derivatives using the Bianchi identities in the 
3 + 1 form [36] 



dtSij =C p E l3 + a DkBipetf - 3£%K. 



d t Bij 



--jOpB^ + 
+KtB t 



kl Y2 rnn 



A7 



(4.16) 



where C denotes Lie derivative, D is induced 3-D covari- 
ant derivative operator, a denotes the lapse, (3 the shift, 
K the extrinsic curvature, and a k = d^hva. The time 
derivative of the metric dt gtj is already known from the 
numerical evolution of the spacetime. Using the above 
equations and applying the chain rule, we compute the 
time derivatives of f and 6: 



(dtgij,d t £ij,d t Bi 



Eqns. (2. 13), (3. 15), (3. 29) 



d t r,d r 



Equipped with all the components of df and dO, we can 
apply Eqs. (3.33) and (3.34) to fix spin-boost degree of 
freedom, finally obtaining the QKT on which we can then 
extract Newman-Penrose scalars via Eqs. (2.4-2.8). 

We note that it may not always be possible to de- 
fine the f and coordinates using Eq. (3.29) for space- 
times with additional symmetries. For example, in ax- 
isymmetric spacetimes admitting a twist-free azimuthal 
Killing vector, ^2 is real, and as a result the 9 coordi- 
nate cannot be computed using Eq. (3.29). In fact, for 
Minkowski spacetimes, we cannot even define the f coor- 
dinate, because ^2 = 0. In such cases, the symmetries of 
the spacetime typically provide a set of preferred coor- 
dinates, which one would naturally adopt in a numerical 
simulation. In our QKT implementation, we presume 
that any such preferred coordinates are adopted, and we 
replace (df,d9) by their simulation-coordinate counter- 
parts when degeneracies occur. 



B. Extrapolation 

We now turn to extracting the asymptotic gravita- 
tional wave content at by using the peeling property, 
i.e., to extrapolation, which necessarily involves informa- 
tion from several spatial slices in the spacetime. Our pro- 
cedure is to shoot a null geodesic affinely parametrized 
by r toward J^ + , monitoring ^4 along the geodesic. The 
best possible polynomial in 1/r is fitted to the result. 
The existence of this polynomial follows from the peel- 
ing property, which is made explicit in Eq. (3.48). Wc 
identify the coefficient of the 1/r term or with the 
radiation content at J? + . 



In contrast to the usual method (extrapolating ^4 
as computed using a tetrad parallel-transported along 
an outgoing null geodesic), note that here we choose 
to extrapolate ^4 (defined using the QKT), which we 
expect to also display the correct peeling behavior [see 
Sec. IIIH3]. In addition, the initial direction of the out- 
going null geodesic is along I, so at the geodesic's starting 
point \E , 4 = if 4, and [Sec. IIIH], at J? + also the outgoing 
null geodesic is along I so that \&4 = \&4. In practice, 
as we integrate along these outgoing null geodesies, we 
monitor the difference between the null vector I tangent 
to the outgoing geodesic and I from the QKT, and we 
find that this difference remains small (cf. Fig. 3 and the 
surrounding discussion). Therefore "J 4 and ^4 are not 
significantly different for the simulations we examined. 
When we extract the ^4 waveform, it converges rapidly 
to its asymptotic value with increasing extraction radius 
[Fig. 16]. 

Selecting the initial tangent of the geodesies to be / 
determines the parameterization of these geodesies upto 
an additive constant B corresponding to the freedom to 
shift the zero point of the affine parameter, r — > r + B. 
The asymptotic waveform is insensitive to the choice of 
the field B. Nevertheless, to provide an exact prescrip- 
tion wc fix B by recalling that in the Kerr limit, the 
affine parameter is just the Boycr-Lindquist f [15]. Wc 
thus choose B on the initial world tube (where we start 
shooting out null geodesies) to be such that r = f. 



C. Sensitivity of QKT method to numerical error 

The numerical implementation of the QKT described 
in this section keeps the computation "as local as pos- 
sible" in the following sense: the bulk of the calculation 
requires only local derivatives and knowledge of the met- 
ric and the extrinsic and intrinsic curvature of the spa- 
tial slice. However, this says nothing about the accuracy 
of our method, which depends on how susceptible our 
method is to numerical noise. 

To begin addressing this issue, we first recall exactly 
how many numerical derivatives are to be taken. Equa- 
tion (3.11), which is used to construct the gravitoelec- 
tromagnetic tensors £ and B, requires i) second spatial 
derivatives of the spatial metric in order to get the in- 
trinsic Ricci curvature, in addition to ii) the first spatial 
derivatives of the extrinsic curvature of the slice. Once 
the gravitoelectromagnetic tensor is obtained and the re- 
sulting curvature invariants / and J are computed, an- 
other derivative is required to compute the gradients of 
the coordinates (f,8) that then fix the Type III free- 
dom of the tetrad. Note that the first step, i.e. the 
computation of the gravitoelectromagnetic tensors only 
requires spatial derivatives, which we can compute spec- 
trally (i.e., inexpensively and accurately, since wc expect 
to observe exponential convergence in spatial derivatives 
with increasing spatial resolution). However, taking the 



gradient of the coordinates constructed out of the cur- 
vature invariants requires both spatial derivatives and 
a time derivative. Fortunately, this time derivative can 
be computed using the Bianchi identities as described in 
Sec. IV A 3, which again reduces the operation to spatial 
differentiation (although here the accuracy of the deriva- 
tives are also limited by the accuracy at which the con- 
straint equations are satisfied). 

What we find in practice is that the higher derivatives 
needed by our QKT method can at places have a signif- 
icantly higher amount of numerical noise than the nu- 
merical derivatives directly used in the actual evolution 
system. This is a significant challenge to our method, 
since SpEC presently evolves the Einstein equations in 
first-order form, i.e., as a set of coupled partial differen- 
tial equations containing only first derivatives in space 
and time. Therefore, the evolution equations themselves 
will only guarantee the existence of one derivative of the 
evolution variables (e.g., of the metric). Constraints show 
convergence which means, among other things, that the 
auxiliary variables (defined during the reduction of sec- 
ond order differential equations to first order) do converge 
to the appropriate metric derivative quantities. How- 
ever, the evolution system, although quite capable at 
constraining the size of numerical error, does not nec- 
essarily force it to be smooth (diffcrcntiablc to higher 
orders) at subdomain boundaries. 

Consider the hypothetical example of adding white 
noise to a smooth analytical metric, such as the Kerr 
metric (2.20). No matter how small the magnitude of the 
noise, it would prevent us from taking derivatives analyt- 
ically. Numerically, under-resolving the high-frequency 
noise would smooth out the data and allow differenti- 
ations to proceed without significantly amplifying the 
added noise; therefore, we expect that filtering (the spec- 
tral equivalent of finite-difference dissipation) would im- 
prove the smoothness of the numerical data and thus 
reduce difficulty in taking higher numerical derivatives. 
However, such filtering can effectively under resolve not 
only noise but also physical information. In other words, 
overly dissipative schemes tend to be less accurate; there- 
fore, the current choice in SpEC is to dissipate as little as 
possible while still maintaining robust numerical stabil- 
ity. This criterion is different from the use of filtering to 
damp out on short time scales any high frequency modes 
that would be produced during an evolution. 

A better approach for reducing non-smooth numeri- 
cal error is to go directly to their source. The lack 
of smoothness in the constraints observed in a typical 
SpEC evolution is partly due to the penalty algorithm, 
which is known to produce convergent but non-smooth 
numerical errors at subdomain boundaries. [See Fig. 4(b) 
for an illustration of the penalty-algorithm induced non- 
smooth error.] Because this non-smoothness converges 
away with increased resolution, our method is observed 
to be viable given a sufficiently high numerical resolu- 
tion; however, it remains to be seen whether "sufficiently 
high" means "significantly higher" than typical resolu- 
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FIG. 4: (a) contours displaying a pulse of high frequency 
junk radiation propagating outward. Ahead of the pulse, the 
geometric coordinate contours are consistent with the Kerr- 
like initial data, but behind the pulse of junk radiation, the 
spacetime settles down to an actual binary inspiral with a 
signature "spiral-staircase" pattern also seen in Fig. 2(b). (b) 
The surface is a 2-D spatial slice containing the symmetry axis 
in a head-on simulation, warped and colored according to f 
value. The sub-domain boundaries are marked out with dense 
black lines, and appear to be a source of non-smooth noise. 
These noisy features are reduced by increasing resolution. 



tions currently in use. Alternatively, improvement to 
non-smooth numerical error could come through the use 
of newer inter-patch boundary algorithms, such as Dis- 
continuous Galerkin methods [37]. There also exists an 
ongoing effort to bring a (currently experimental) first- 
order-in-time, second-order-in-space version of SpEC [38] 
into a state suitable for accurate gravitational-wave pro- 
duction, with the hope of added efficiency and of achiev- 
ing numerical error of higher differential order. Such pos- 
sibilities as these, however, are future work, well outside 
the scope of this paper. 

Lastly, we consider the non-smooth noise sensitivity of 
our QKT quantities from another point of view: it can 
be used as a diagnostic of high-frequency, non-smooth 
numerical error. For instance, one source of non-smooth 
constraint violation in numerical simulations is the high- 
frequency, spurious "junk" radiation present at the be- 
ginning of numerical simulations (because of how the ini- 
tial data are constructed), which poses a particularly dif- 
ficult numerical problem. The frequency of these modes 
is of O(M), orders of magnitude higher than that of the 
orbital motion (and the associated gravitational waves). 
This makes resolving the junk radiation a difficult task. 
In the effort to reduce junk radiation, the geometric coor- 
dinates can be used as a visualization tool. Fig. 4(a) is an 
illustration of how the 9 contours, plotted as a function 
of code coordinates, react to the junk travelling through 
the grid, while adjusting themselves to reflect a more re- 
alistic spacetime. By comparing the difference ahead and 
behind the easily identifiable junk pulse in Fig. 4(a), one 
gets a glimpse of the missing pieces in the initial data. 
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FIG. 5: A Kerr-Schild black hole with J/M 2 = 0.5, with the 
coordinates translated a distance 9M along x axis. The red 
sphere indicates location of the black hole's horizon, and the 
black cross indicates the coordinate origin. The colored circles 
are constant geometric radius f contours; these demonstrate 
the ability of our geometric coordinates to select an origin 
based on the Coulomb potential of the QKF, i.e., an origin 
which reflects the gravitational curvature of the spacetime. 
Also shown are the spatial projections of the n direction of 
the coordinate and quasi-Kinnersley tetrads, at points on a 
narrow strip marked by a grey ring. The black lines indicate 
the n direction associated with the coordinate tetrad, which 
point toward the coordinate origin, and the red lines/arrows 
are QKT n directions that identify the black hole as geometric 
origin, away from which the gravitational waves travel. 



NUMERICAL TESTS OF THE QKT 
SCHEME 



We now consider several numerical tests used to gauge 
the effectiveness of our proposed QKT scheme for wave- 
form extraction. Most of these tests are motivated by 
analytic solutions and are used to verify that our choices 
of geometric coordinates and the QKT are yielding the 
expected results. These tests broadly fall into two classes: 
i) non-radiative spacetime tests and ii) radiative space- 
time tests. Each will be considered in turn in the follow- 
ing subsections. 



A. Non-radiative spacetimes 

1. Kerr black hole in translated coordinates 

The spacetime in this test is a Kerr black hole in Kerr- 
Schild coordinates, but the coordinate origin is trans- 
lated away from the black hole along x or z axis. Here 
we work in units of the black hole mass, and the dimen- 
sionless spin is J/M 2 = 0.5 pointing in the z direction. 
Tetrads determined only by our simulation coordinates 
[see Eqs. (4.2)-(4.8)] would not be aware of the transla- 
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FIG. 6: (a): The L2 norm [Eq. (5.1)] of Newman-Penrose 
scalar ^4 computed (between radii 50M and 140M) for the 
Kerr-Schild black hole in translated coordinates on the coor- 
dinate tetrad, (b): The Newman-Penrose scalar ^4 computed 
in the quasi-Kinnersley tetrad. Note the L2 norm is eleven 
orders of magnitudes smaller for the QKT when compared to 
that of the coordinate tetrad, showing that the QKT correctly 
adapts to the underlying curvature of the spacetime. 



tion, and the spatial projection of n would point toward 
the coordinate center instead of the black hole itself. In 
contrast, the QKT should adjust to the displaced origin, 
picking up the true geometrical origin of the gravitat- 
ing system determined by the Coulomb potential of the 
QKF. Figure 5 shows the direction of spatial projection of 
n and n associated with the two tetrads. The QKT iden- 
tifies the black hole at the center of the circular shape, 
as do the geometrically motivated coordinate f. 

Figure 6 compares ^4 extracted using the coordinate 
and quasi-Kinnersley tetrads, respectively, using the so- 
called "L2 norm" as a measure. The L2 norm of a quan- 
tity X is defined here as 



\ 



N u 

E 



Ntnt 



(5.1) 



where Xi are the spectral collocation points of a pseudo- 
spectral grid and N to t is the total number of points. The 
present study uses four spherical shells between radii 50M 
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and 140M with N tot w 4 x 45 3 collocation points. The 
QKT correctly produces vanishing (up to numerical 
round-off error), while the coordinate tetrad fails to iden- 
tify the correct out-going direction and as a result misin- 
terprets ^2 as gravitational radiation content in ^4. (We 
observe similar behavior for ^o-) Using such a coordinate 
tetrad in a simulation with a displaced center will result 
in spurious effects being picked up in the extracted ra- 
diation, of a magnitude not necessarily smaller than the 
physical gravitational wave content of the spacetime. 

In a simulation of a dynamical spacetime, a similar ef- 
fect should be expected when the "center of mass" (e.g. 
in a Newtonian approximation) of the system does not 
coincide with coordinate center. For example, consider a 
binary merger of unequal mass holes with the coordinate 
origin placed at the midpoint between the black holes; 
^4 extracted at finite radii would pick up a slowly vary- 
ing offset at an integer multiple of the orbital frequency, 
and this contribution would complicate the extrapolated 
waveform. 



2. A Schwarzschild black hole with translated coordinates 
and a gauge wave 

We further explore the effects of coordinate choice or 
gauge by introducing a time dependent gauge wave into a 
Schwarzschild solution whose origin has been translated 
by a constant amount. The resulting metric components 
now have an explicit time dependence, and we expect the 
coordinate tetrad to produce a false gravitational wave 
signal, even though the Schwarzschild spacetime is static 
and emits no physical radiation. 

The exact analytic solution we use for this test is 
constructed from the Schwarzschild solution in ingoing 
Eddington-Finkelstein coordinates. We then apply a 
time-dependent coordinate transformation that yields a 
metric of the form 



ds 2 = - (H 
+ 2(1 

+ (H 
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C) 
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C 



dtdr 



2M 
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C 



dr 2 + r 2 dfl 2 



(5.2) 



where C (r, t) is the radial waveform of the introduced 
gauge wave. For our test we select generically chosen 
parameters 



C = 0.7sin(0.03(t + r) +3.1), M = \ 



(5.3) 



Note that again we translate the black hole off the coor- 
dinate origin by a constant amount (here r = 20M) as 
described in the previous subsubsection. 

Figure 7 shows spin-weighted spherical harmonic ex- 
pansion coefficients \J/^' m - ) of ^4, computed using the co- 
ordinate and quasi-Kinnersley tetrads. While only the 
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FIG. 7: Spherical harmonic coefficients of Newman-Penrose 
V?4 computed for the gauge wave solution, on a sphere of 
radius r = 120M centered on the black hole, (a): ^4 extracted 
on the coordinate tetrad, (b): <&4 extracted on the QKT. 
Only the three largest (I, m) spin- weighted spherical harmonic 
modes (up through / = 35) are shown. Note the scaling on 
the two figures differ by nine orders of magnitude. 



three largest amplitudes are shown, we have computed all 
amplitudes up through I = 35. These scalars are com- 
puted on a sphere at a radius of 120M from the black 
hole, with the poles of harmonics aligned with the direc- 
tion in which the black hole is shifted. As expected, the 
waveform extracted using the coordinate tetrad picks up 
a time dependence associated with the gauge wave, while 
the QKT returns vanishing values, correctly identifying 
the static spacetime solution. 

In the generalized harmonic form of the Einstein field 
equations, the gauge may be set by the covariant wave 
equations 



Ux a = H a 



(5.4) 



where H is either a specified or evolved source func- 
tion [39-42] . It is thus probable that gauge modes similar 
to the one considered in this example may be present in 
fully dynamical simulations. Consider a gauge wave that 
generates a deviation between the coordinate tetrad ba- 
sis vectors {I, n, m, rr%] and their counterparts in the 
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QKT. Such differences can be represented by a sequence 
of type II, I and then III transformations parameter- 
ized by the time dependent transformation parameters 
b(t), a(t) and A(t) that appear in Eqs. (2.15), (2.14) and 
(2.16) respectively. If we restrict ourselves to asymptotic 
regions where ^4 dominates over other NP scalars, then 
according to Eqs. (2.18), (2.17) and (2.19), we have that 
to leading order in a and b the coordinate "J 4 is given by 



* 4 = {l + Aa{t)b{t))A- 2 {t)^i. 



(5.5) 



If the gauge wave falls off when we move away from the 
source region, then we may have 



(l + 4a{t)b(t))A- 2 (t) -> 1 



(5.6) 



and its effect can in principle be extrapolated away. How- 
ever, for some cases, such as a plane gauge wave, the time 
dependent perturbation introduced into ^4 could persist 
in the extrapolated waveform. Therefore, minimizing any 
such gauge-dependent content in ^4 extracted at finite 
radii is preferable to relying on extrapolation to remove 
them; some pathological gauge modes might not fall off 
sufficiently quickly with radius. 



B. Radiative spacetimes 

Having observed that the QKT correctly reflects the 
curvature content of non-radiative spacetimes, including 
in the presence of a gauge wave, we next apply the QKT 
to spacetimes emitting gravitational radiation. In this 
subsection, we verify that the scheme is consistent with 
analytic perturbation theory results. 

The QKT by construction reduces to the Kinnersley 
tetrad in the Kerr limit. Therefore, if we perturb a 
Kerr black hole by a small amount, the ^4 computed 
on the QKT should reproduce the analytic perturbation 
theory results computed on the Kinnersley tetrad asso- 
ciated with the unperturbed Kerr background. Verify- 
ing this correspondence provides us with the means to 
quantitatively test whether the QKT extracts the correct 
waveform and that we have all normalization conventions 
implemented correctly. The idea of ensuring the corre- 
spondence between the computed waveform and the per- 
turbation theory results is what motivated the authors 
of Rcf. [9] to adopt transverse tetrads in the first place; 
Chandrasckhar [15] also used the transverse tetrad in his 
metric reconstruction program, where he explicitly com- 
puted the perturbed tetrad and curvature perturbations 
on the tetrad, obtaining the expected correspondence. 
For simplicity, here we perturb a Schwarzschild black hole 
with an odd-parity Rcggc-Wheeler-Zerilli (RWZ) pertur- 
bation, as described in Ref. [43]. 

We start with a background Schwarzschild metric in 
Schwarzschild coordinates expressed in the standard form 
of [43] 




(b) 



FIG. 8: Sft(\p4) resulting from a traveling-wave perturbation 
on the equatorial plane of the computational domain. Panels 
(a) and (b) correspond to results obtained with and without 
the coordinate transformation (5.20). We use ^4 , ^4 and 
\l/4 to respectively denote the analytical result and the val- 
ues computed on the coordinate and quasi-Kinnersley tetrads. 
The height of the surface in the vertical direction indicates the 
value of K(^&4), the solid grey surface denotes Sft(*& 4 ), the red 
wire-frame Sft(\p4) and the black dots ^(^4). The amplitude 
of the red wireframe has been suppressed by a factor of 10 3 
in (a) so that it fits into the figure. The suppression factor 
has not been applied to Panel (b). 
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We then introduce I = 2, m = ±2 radiative perturba- 
tions. The full RWZ formalism giving the explicit calcu- 
lation of the perturbed metric is expounded concisely in 
Appendix A of [43] ; it turns out that the construction of 
the perturbed metric and the associated perturbed cur- 
vature quantities (such as iff 4) hinges on one function, 
the RWZ function Z, which obeys the RWZ equation 
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In the RWZ equation, the coefficients Cj are functions 
of a, j3 and 7, which for our chosen values [Eq. (5.8)] 
become 
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(-a 2 + )dt z + 2 1 2 Pdtdx + ~i 2 dx 2 



Given Z, the analytic solution ^4 for the gravitational 
wave content in the spacetime can be computed to first 
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FIG. 9: Testing the QKT's ability to recover perturbation theory results. An I = 2, m = ±2 perturbation, with magnitude 
Ac, is added as a function of retarded time in addition to an I = 1, m = perturbation, and finally a cubic coordinate 
distortion is applied as described in Eq. (5.20). The spherical harmonic coefficients of the coordinate-tetrad ^4 and QKT $4 
are then extracted at Boyer-Lindquist radius 95M and compared with the analytic results *P 4 . (a): The magnitude of the ^4, 
(Z,m) = (2,0) mode as a function of spin parameter a, with the radiative (2, ±2) perturbation held fixed at A c = 10 _9 M. (b): 
Exploring the effect of cubic rescaling. The rescaled first component of f vector, f , is plotted against the original component 
r . The identity map (f 1 = r 1 ) is given for comparison, (c) and (d): The real and imaginary parts of the -2V22 coefficient 
vs the amplitude A c of the radiative perturbation. The I = l,m = mode amplitude is chosen so that the resulting angular 
momentum perturbation is held constant at a — 0.001M. 



order using 
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(5.14) 



where Ci = ^J{1 - 1)1(1 + I) (I + 2)/4, the operator A is 
n a V a with n being a null direction associated with the 
background Kinnersley tetrad, and 7 and are spin co- 
efficients associated with the same background tetrad, 
which for our case are [15] 
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A 



r-M 
2S 



(5.15) 



where A, E and p are defined in Eq. (2.21). 

Note that in the discussion that follows, ^4 denotes 
the analytic result while ^4 and ^4 are, respectively, the 
computed values on the coordinate and quasi-Kinnerslcy 
tetrads in the numerical implementation. 



In order to solve the linear second order partial dif- 
ferential equation (5.9) for Z, an initial value and time 
derivative for Z must be specified. For our investigation, 
we make use of a traveling- wave perturbation of the form 

87 



dt 



(5.16) 



where r* is the usual tortoise coordinate defined by 
dr^/dr = r/(r - 2M), while t = 0, u = 0.1 and A c 
is a constant initial amplitude. For our test, we also set 
M = 1. This perturbation is graphically depicted in Fig. 
8. The waveform constructed from the perturbation has 
the classical profile for ^4, often observed during numer- 
ical binary black hole mergers; this is to be expected, 
since I = 2 is the dominant mode contributing to the 
gravitational radiation emitted by a binary. 

Next, we numerically compare the coordinate-tetrad 
^4 and the QKT ^4 with the analytic perturbation- 
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theory result In this test, we adopt the Boyer- 

Lindquist coordinates; therefore, the corresponding coor- 
dinate orthonormal tetrad [see Eq. (4.2)] happens to co- 
incide with a Kinnersley frame of the background space- 
time [although it is boosted with respect to the Kinner- 
sley tetrad in Eqs. (2.22)- (2.24)]. 

To illustrate this more clearly, observe that the stan- 
dard coordinate tetrad we constructed in Sec. IV A 1 re- 
sults in orthonormal vectors T a and N a that are respec- 
tively 



rpa 



1 



-, 0, 0, 



N a = [0, a, 0, 0] (5.17) 



when expressed on the coordinate basis, where a = 
y/T— 2M Jr is the lapse of the background metric. The 
resulting null vector l a constructed according to Eq. (2.1) 
is 



l a = -=(T a + N a ) = -4= 



-,a,0,0 
a 



(5.18) 



In the static limit, the Kinnersley tetrad [Eqs. (2.22) and 
(2.23)] reduces to 



r = 



\, i, o, o 

cr 



, n a = 1 [1, -a 2 , 0, 0] (5.19) 



so l a = \/2l a /a and there exists a relative boost factor 
of A = aj\pl between the coordinate and Kinnersley 
tetrads. 

Therefore, to account for the difference, we will multi- 
ply the coordinate-tetrad ^4 by (1 — 2M/r)/2 throughout 
this subsection to facilitate comparison with analytical 
and QKT values. With this adjustment, the extracted 
quantities ^4 and V&4 both match the analytically cal- 
culated j£ 4 from perturbation theory. These results are 
graphically depicted in Fig. 8(b). 

Next, we explore the gauge dependence of the QKT 
result. To this end, we introduce a coordinate trans- 
formation into some other gauge. As a result, the co- 
ordinate tetrad associated with the new "non-privileged 
gauge" differs from the QKT, resulting in a mismatch be- 
tween the coordinate "J 4 and the analytic perturbation 
theory result £4 [see Fig. 8(a)]. The QKT ^4 imple- 
mented in the code should then be able to recover the 
analytic result. As an illustrative example, for a gauge 
transformation we choose a cubic rcscaling of the spatial 
coordinates, which takes the radial vector r a expressed 
on a Cartesian coordinate basis defined in Sec. IV A to a 
vector with components f a using the equation 



using both the coordinate tetrad and the QKT, com- 
paring the results with the analytical \&4 calculated in 
the Boyer-Lindquist coordinates and visually portrayed 
in Fig. 8. When no coordinate distortions have been in- 
troduced [Fig. 8 (b)], the coordinate tetrad and the QKT 
both generate 'J 4 that matches the analytical prediction. 
When we apply the cubic coordinate distortion described 
above however, the coordinate tetrad result deviates from 
_^4, but the QKT still recovers the analytical value [Fig. 8 
(a)]. 

In addition to the coordinate transformation we also 
add a I = 1, m = mode [explicit expressions for metric 
perturbation due to this mode can be found in Appendix 
Ala of Ref. [43]]. This mode should make no contribution 
to the detected radiation in ^4, but should introduce a 
small angular momentum perturbation affecting the spin 
of the spacetime, thus avoid degeneracy in 9. The am- 
plitude of the perturbation is usually set so that the spin 
of the resulting spacetime is a = J/M = 0.001. After 
imposing the coordinate distortion and adding the I = 1 
m = mode, we extract -■^Xim coefficients of ^4 and 
^4 on a sphere of Boyer-Lindquist radius 95M from the 
black hole. We begin by exploring the effect of the 1 = 1, 
in = mode. Figure 9 (a) shows the effect of increasing 
the strength of the I = 1, m = perturbation on the 
I = 2, m = mode of extracted waveform; recall that we 
did not introduce an I = 2, m = mode into the met- 
ric perturbation. The coordinate quantity "J^ 2,0 ^ shows a 
constant value possibly originating from the cubic coordi- 
nate distortion. The quantity ^^' \ on the other hand, 
shows a strong dependence on the I = 1 perturbation am- 
plitude; this could be due to the fact that the projection 
onto spherical harmonics, as opposed to spheroidal har- 
monics, is no longer correct when the perturbing spin is 
introduced. The effect is however small when compared 
to the magnitude of the I = 2, m = 2 modes. 

The lower panels of Fig. 9 explore the effect of increas- 
ing the I — 2, m = ±2 perturbation amplitude A c on the 
extracted I = 2, m = ±2 modes. In general the QKT 
shows very good agreement with analytical result over a 
range of perturbation magnitudes as desired, while the 
quantities computed on the coordinate tetrad disagree 
significantly with the analytic perturbative result. 



VI. APPLICATION OF THE QKT TO 
NUMERICAL SIMULATIONS OF BINARY 
BLACK HOLES 



R 2 



\r-r \* (r a -r a )+r a . (5.20) 



where we choose rg = (5M, 0, 0), R = 100M, v = 1 and 
v = 1.1. Panel (b) of Fig. 9 compares the coordinates 
before and after rcscaling. 

We now calculate the perturbed metric in the new dis- 
torted coordinates and extract the gravitational waves 



We now turn to exploring the properties and effective- 
ness of the QKT scheme when applied to more generic 
numerical simulations involving the collision of two black 
holes. We consider two examples: a circular inspiral 
of two equal-mass, nonspinning black holes [Sec. VI A] 
and the head-on collision of two nonspinning, equal-mass 
black holes [Sec. VI B]. 
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A. Equal-mass, nonspinning binary-black-hole 
inspiral 

In this subsection, we apply our QKT method to a 
fully dynamical simulation of two equal-mass, nonspin- 
ning black holes that inspiral through 16 orbits, merge, 
and ring down. We summarize some of the physical pa- 
rameters of this simulation in Table I (which is a repro- 
duction of Table II of Ref. [44]); further details of this 
simulation and the numerical method used are given in 
Ref. [441 and the references therein. 



Initial orbital eccentricity: 
Initial spin of each hole: 
Duration of evolution: 
Final black hole mass: 
Final spin: 



e ~ 5 x 1(T J 
Si/M 2 < 10" 7 
AT/M = 4330 
M f /M = 0.95162 ± 0.00002 
Sf/Mf = 0.68646 ± 0.00004 



TABLE I: Physical properties of the equal-mass, nonspinning 
binary-black-hole inspiral reported in Ref. [44] . Here M is the 
sum of the Christodoulou masses of the initial holes, and M/ 
is the Christodoulou mass of the final hole. 

We examine two aspects of the QKT that we have con- 
sidered in previous sections: i) that the direction I iden- 
tified by the QKT corresponds to the wave-propagation 
direction (as discussed in Sec. Ill A) and the implica- 
tions this has for the geometric coordinates f and 9, and 
ii) that the falloff rates of the Newman-Penrose scalars 
are consistent with the peeling property (as described in 
Sec. IIIH). 







FIG. 10: Geometrical coordinates obtained from the QKF \&2- 
(a) Contours of r on a slice of the null preferred characteris- 
tic surface (PCS) generated by the geodesic developments of 
I. (b) Contours of f on a constant-simulation-time Cauchy 
slice, (c) Contours of on the same null surface as (a), (d) 
contours in the same Cauchy slice as (b). 



mass quadrupolc approximately given by the Newtonian 
relation 



1. Wave-propagation direction 

If (as claimed in Sec. Ill A) the vector I associated 
with the QKT correctly identifies the out-going wave- 
propagation direction, one would expect that as one fol- 
lows the wavefront out to infinity, the spacctime cur- 
vature along this trajectory and the associated derived 
quantities should become quite simple. To illustrate 
this, we consider an S 2 coordinate sphere in the origi- 
nal Cauchy slice and identify the correct null direction 
I associated with the QKT at each point. We then in- 
tegrate the null geodesic equations outward to produce 
a null hyper-surface and consider this null hyper-surface 
to be a new slicing and a preferred characteristic sur- 
face (PCS) of the spacetime. Identifying the geometri- 
cal coordinates f and 6 within this slicing, we plot their 
contours in the left-hand panels [plots (a) and (c) respec- 
tively] of Fig. 10. For comparison, we also show f and 
6 computed within the original Cauchy slice. Note the 
simplicity of the computed geometric quantities within 
the PCS associated with the out-going wavefront as op- 
posed to the corresponding quantities computed within 
the Cauchy surface. 

The structure observed within the Cauchy surface can 
be understood as follows. The holes generate a rotating 
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where M is total mass of the binary, R is the separation 
between the two black holes, q is the location of one of 
the black holes, and the choice of coordinates is such 
that t he other black hole is located at —q. In the matrix, 
O = y/ Mj R 3 is the orbital angular velocity and t' = t—r. 
This quadrupolc moment deforms the r contour into an 
ellipsoid (or peanut shape when closer to the two holes 
[see Fig. 2 (a)]), while its time dependence causes the 
orientation of the ellipsoid to rotate at a frequency of 
2il. 

On the PCS, the structure is much simpler. The in- 
ner contour sets the basic shape for constant f surfaces, 
which are roughly ellipsoidal. These surfaces then ex- 
pand, retaining their orientation as the distortion is prop- 
agated outward at the speed of light along the wavefront. 
Figure 10 (a) shows a concentric pattern of f contours 
on the null hyper-surfaces, in contrast to the rotating 
contours on a spatial Cauchy hyper-surface that slices 
through many PCSs, as depicted in Figure 10 (b). The 
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Recall that in fixing the spin-boost or Typo III freedom 
of the QKF to obtain the QKT, the I vector was scaled 
so that ^2 oc (r)~ 3 . The very sharply defined peak at 
—3 in Fig. 11 provides direct numerical evidence that the 
relation r oc f [cf. Sees. IIIH and IV B] remains valid at 
leading order for the considered range of the computa- 
tional domain. 

Figure 11 also indicates that |^o| and |^4| scale as r -5 
and r _1 , respectively, as expected from Eq. (3.48). Here, 
the peaks are not as sharply defined, since we do not by 
construction enforce the power-law scalings of an d $4 
(as we do for $2)- 



B. Head-on nonspinning binary merger 



FIG. 11: The distribution of power-law falloff rates of 
Newman-Penrose scalars against affine parameter r. The 
three concentrations (colored blue, green and red) from left 
to right indicate falloff rates of Wo, ^2 and ^4, respectively. 
The vertical axis indicates the number of geodesies (totaling 
3510) with falloff rate falling inside bins of width 0.08, and 
the number of geodesies for the center-most bins are shown. 



angular 9 coordinates similarly display a relative simplic- 
ity on the PCS, taking on the shape of a slightly deformed 
(squashed sideways) cone. Figures 10 (c) and (d) show 
the 9 surfaces on a PCS and in a spatial slicing respec- 
tively; the orientation of the deformed constant 9 cones 
is independent of the distance to black hole in the PCS, 
but rotates around when moving outwards on the spatial 
slice, forming a "spiral-staircase" pattern. 



2. Peeling property 

Next, we explore the falloff rate of the Newman- 
Penrose scalars computed on the QKT as one moves 
outward along the PCS generators (i.e., along the null 
geodesies tangent to the QKT I where they originate). 
This rate allows us to quantify to what extent the QKT 
obeys the peeling property derived in Sec. IIIH and to 
what extent the computed quantities are suitable for use 
in the extrapolation procedure prescribed in Sec. IV B. 

To this end, we start with 3510 null geodesies from the 
grid points of a mesh (of 351 points) covering a sphere 
of radius f ~ 150M surrounding the source region. Over 
a small time interval of 10M, a new set of geodesies are 
shot off every 1M. The affine parameter r is initially set 
to f and the geodesies are evolved for around 150M. The 
Newman-Penrose ^i's are recorded at intervals of At = 
1M along the geodesies. A histogram of the best fits for 
the power-law falloff (i.e., of the slopes of the ln(|^i|) vs 
ln(r) graphs) for the 3510 geodesies are plotted in Fig. 
11. 



To further examine the properties of the QKT and the 
geometrical coordinates we now take a detailed look at 
the numerical simulation of a head-on merger. The phys- 
ical parameters of the simulation are given in Table II. 



Initial separation: 
Initial spin of each hole: 
Duration of evolution: 
Final black hole mass: 
Final black hole spin: 



d/M = 20 
S/M 2 < 2 x 10" 12 
AT/M = 600 
Mf/M = 0.987 ± 2 x 10 



Sf/Mf 



3. x 10"' ± 2 x 10 



-3 
7 



TABLE II: Physical parameters of the head-on binary-black- 
hole merger considered in Sec. VI B. Here M is the sum of the 
initial black hole Christodoulou masses, and all initial quan- 
tities are measured at the initial time t = of the simulation. 



The axisymmetric head-on collision of two nonspinning 
black holes has been studied extensively [45-49] ; in many 
respects, these collisions serve as a simple, strongly non- 
linear test of numerical relativity codes. The existence 
of a twist-free azimuthal Killing vector on this spacetime 
implies that the metric does not explicitly depend on 
the azimuthal coordinate <p defined about the symmetry 
axis and that the angular momentum of the spacetime 
is zero. We note that because of the symmetry of this 
configuration, the Coulomb potential associated with the 
transverse frame is real and thus that only one geometric 
coordinate, the radial coordinate f, can be determined 
from it. Therefore, we fix the latitudinal coordinate to 
the simulation coordinate 9. 



1. Geometric radial coordinate 

We now explore some of the properties of radial coordi- 
nate f, the emitted radiation profile, and the waveform. 
We show contour plots of the f coordinate at various 
times near merger in Fig. 12. The characteristic peanut 
shape expected from the merger event is clearly visible, 
and surfaces of constant f coordinate trace both the in- 
dividual apparent horizon surfaces at early times and the 
final apparent horizon surface at late times. Far from the 
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FIG. 12: Here we show the evolution ol f contours during 
the merger in the near zone; the black surfaces in these plots 
are the apparent horizons, and at t = 112M, the common 
apparent horizon forms. 

source, constant f surfaces become roughly spherical, in- 
dicating that there the geometrical concept of radius and 
the gauge choice for the radial coordinate in the simula- 
tion coincide well. For the head-on collision, as well as in 
the more dynamical spacetimes depicted in Fig. 2, plot- 
ting surfaces of constant f turns out to be a useful tool 
for visualizing the spacetime geometry in a way that cor- 
responds to an intuitive feel of the Coulomb potential's 
behavior. 



2. Gravitational waveform 

For twist-free axisymmetric spacetimes, one can show 
in general [45] that if the imaginary part of the tetrad 
null vector m or E 3 (as defined in Sec. II A) has the 
same direction as the azimuthal Killing vector, then ^4 
expressed on this tetrad is real. Figure 13 depicts ^(^4) 
for the head-on collision presently under consideration. 
Note the absence of radiation along the symmetric axis 
in Fig 13; this is a feature we will examine further later 
in this section in the context of PNDs. 

We show two possible spherical harmonic decomposi- 
tions of ^4 in Fig. 14. Panel (a) corresponds to the case 
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FIG. 13: A snapshot (at t = 242. 25M) of latitudinal distri- 
bution of radiation emitted by the head-on collision. The col- 
oring is according to |5R(*I/4)|, with large values corresponding 
to darker color. The disk is a vertical slice of the computa- 
tional domain with the thick red line denoting the symmetry 
axis. 



where the azimuthal Killing vector determines the 9 = 
direction; because the symmetry-axis corresponds to the 
9 = direction, axisymmetry implies that there are no 
I = 2, m = ±2 modes in the spherical harmonic decom- 
position, only ttl — modes exist, and of those the 1 = 2, 
m = mode makes the dominant contribution. 

On the other hand, if one relabels the 9 and (f> coor- 
dinates on the extraction sphere, the spherical harmonic 
decomposition of the same waveform is very different. 
Panel (b) of Fig. 14 instead chooses the 9 = line to 
be orthogonal to the axisymmetry axis rather than along 
it [as in panel (a)]; a significant I = 2, m = 2 mode 
appears. This is a simple example illustrating the well- 
known fact that unless a clear prescription for the pre- 
ferred axis of a simulation is given, the I = 2, m = 2 mode 
is an ambiguous description of the radiation. Solutions 
to this problem for generic black hole binary simulations 
(which include precession) have been proposed in liter- 
ature. For example, one may first choose a "radiation 
axis" [50, 51] to maximize the component of the angular 
momentum along itself and secondly choose a preferred 
rotation about that axis [52]. Although not yet fully 
explored, our geometric coordinates suggest an alterna- 
tive resolution. Namely we can use the extrema of the 
computed geometric 9 coordinates to identify the polar 
regions of a simulated spacetime. 

Another question especially relevant to wave extrac- 
tion is how rapidly the waveform computed from \p4 
in the computational domain converges to "the" cor- 
rect asymptotic waveform. In Sec. Ill H, we argue that 
asymptotically the QKT quantities on the correct out- 
going geodesies (as described in Sec. IV B) should con- 
verge very rapidly to the desired result. We now explore 
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FIG. 14: (a) The <E , 4 waveform for the head-on simulation ex- 
tracted at coordinate radii of r = lOOAf and r = 150M. For 
the waveform at r = 100M, all -lYim modes up to I = 35 
are shown. Only K^^ 2 ' ^) and K^^ 4,0 ') are discernibly 
non- vanishing, with the former clearly dominating. For the 
r = 150Af waveform, only K^^ 2,0 ') is shown, it is shifted 
temporally and rescaled by ~ 1.5 so that its maximum peak 
location and magnitude match those of the r — 100M wave- 
form. The difference between the two waveforms is shown 
in the top-right inset, (b) The spherical harmonic decompo- 
sition for ^4 at r = 100M obtained by choosing the poles 
of harmonics along a direction orthogonal to the symmetry 
axis. Multiple modes are visible, and a few with the largest 
magnitudes are labelled. 



this statement quantitatively for the emitted radiation 
on the equatorial plane in the head-on binary-black-holc 
merger we are considering. The goal is to determine at 
which radius a reliable approximation of the asymptotic 
waveform is attained. 

To locate a good cut-off radius for wave extraction, 
consider a set of non-inertial observers hovering at dif- 
ferent fixed spatial (simulation) radii in the equatorial 
plane. For each observer, in Fig. 15 we record the f\E , 4 
value as a function of time and plot the resulting curves, 
with the origin shifted so that the central maxima of all 
the curves coincide at around t = 140M. For clarity, we 
have divided the curves into two sets, those originating at 
r < 30M and at r > 30M, which we display in panels (a) 



and (b) , respectively. In both panels the curve traced out 
at r = 160M is given for comparison, and we will refer 
to it as the reference waveform. Panels (c) and (d) show 
the absolute difference between the curves extracted at 
the various interior points and the reference waveform. 

In Fig. 16, we plot the fractional difference between 
reference waveform and the interior waveform as a func- 
tion of extraction radius. The quantity plotted is the Li 
norm of the absolute difference between the two wave- 
forms divided by the Li norm of the reference waveform. 
The L 2 norm is defined here as 



L2[f] 



p{t)dt. 



(6.2) 



/*=125 



Comparing Fig. 16 and panel (b) in Fig. 15, we find that 
for radii greater than r = 40Af the extracted waveform 
corresponds closely, within ~ 5%, to the reference wave- 
form. Figure 16 quantifies this further: for r < 40M, the 
errors in the extracted \&4 waveform are large but the con- 
vergence to the reference waveform is super-exponential, 
while for r > 40M, the errors converge exponentially. 
This provides quantitative justification for the rapid- 
convergence claims we made in Sees. Ill H with respect to 
the Newman-Penrose scalars calculated on the QKT. We 
conclude that the radius f = 40M appears to be a good 
minimal extraction radius for QKT quantities for head- 
on binary-black-hole collisions. It would be interesting 
to explore, using a similar analysis, whether the expo- 
nential convergence properties and the value for a good 
minimal extraction radius change considerably when ap- 
plied to generic spacetimes (i.e., to spacetimes with less 
symmetry than the head-on collision we consider here). 



3. Principal null directions 

We conclude our investigation of the spacetime associ- 
ated with the head-on collision of two black holes by ex- 
ploring the behavior of the PNDs on and off the axis. As 
we have observed and shown graphically in Fig. 13, no ra- 
diation is emitted along the symmetry axis. This lack of 
radiation suggests that the PNDs do not all converge into 
a Type N pure radiation configuration as seen in Fig. 3 
(c); instead, we would expect the PNDs to remain in a 
Type-D configuration, with two pairs of PNDs trapped at 
antipodal points of the ant i- celestial sphere. On the axis, 
the only non-vanishing Newman-Penrose scalar on the 
QKT is The four solutions in Eq. (3.54) then divide 
into a pair whose value diverge and a pair that approach 
zero. The diverging solutions give us two PNDs pointing 
along the n direction, while the vanishing solutions coin- 
cide with I. This situation is depicted in Fig. 17 (a): the 
two PNDs at the bottom of the sphere point away from 
I (represented by the radial line) and do not converge 
onto the other two PNDs that coincide with I at the top 
of the anti-celestial sphere. All the PNDs on the top of 
the sphere form angles smaller than 4 x 10 _7 7r with the 
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FIG. 15: Panels (a) and (b): The gravitational-wave signal ^(^4) curves extracted at several fixed spatial (simulation) 
coordinates. Panels (c) and (d): The absolute difference between these waveforms with a reference waveform computed at 
r = 160M. 
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FIG. 16: Exponential convergence of the QKT waveform ex- 
tracted in the interior to the reference waveform measured at 
r = 160M. This plot shows the L2 norm [defined in Eq. 6.2] of 
the absolute difference between the waveform at a particular 
extraction radius and the reference waveform at r = 160A/ 
normalized by dividing by the L2 norm of the reference wave- 
form. Also shown is the red fitted exponential curve for radii 
r > 40M, which has a slope of —0.011. 



spatial projection of Z, while the other two PNDs form an- 
gles with the spatial projection of n that are smaller than 



4 x 10~ 6 7r. For comparison, in Fig. 17(b), we show the 
PND behavior on the celestial sphere as one moves along 
a geodesic that points away from the symmetry axis. The 
geodesic shown in this plot starts out at an orientation of 
0.5087T to the symmetry axis. (A similar plot is made in 
Fig. 3(c) where we showed a geodesic starting at 0.3967T 
to the symmetry axis). In these two cases, the presence 
of radiation causes the two PNDs pointing away from I 
to converge onto the other two PNDs surrounding I as 
one moves outward along the geodesic. The dominant 
rate of convergence is 1/t [cf. Eq. (3.53)]. 

The existence of "critical" directions as demonstrated 
here in the special case of axisymmetric spacetimes, is a 
generic feature of all dynamical spacetimes: it is a topo- 
logical necessity [see [53] and page 173 of [34]]. Specif- 
ically, in Ref. [53], the authors explain this feature as 
follows: In the asymptotic region, gravitational radia- 
tion is transverse and can be represented by tendex and 
vortex lines tangent to spheres of constant f. Then, the 
Poincarc-Hopf theorem dictates that there must be loca- 
tions on the sphere where the tendicity associated with 
the two transverse eigenbranches of the gravitoelectric 
tensor £ become degenerate. The trace-free property of 
£ then further constrains the tendicity to be zero — i.e., 
requires the gravitational radiation to vanish at the crit- 
ical points. 
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FIG. 17: A graphical representation of the peeling-off be- 
havior of the PNDs along out-going null geodesies for the 
head-on collision. The null geodesies start at r — 30M at 
time t = 152M. The tangent £ to the geodesies coincide with 
I along the geodesic and is denoted in the figures by black 
radial lines. The geodesic shown in panel (a) travels out- 
ward along the symmetry axis where no gravitational radia- 
tion is emitted. The geodesic shown in panel (b) starts off 
on the equatorial plane, and all PNDs converge toward the 
out-going I direction pointing toward the front of the sphere. 
In this figure, the arrows indicate the movements of the PNDs 
as one travels further along the geodesic. Darker coloring in- 
dicates points calculated at larger affine parameter along the 
geodesic. 

Another useful characterization of dynamical numer- 
ical simulations is a measure of how rapidly the space- 
time settles down to Petrov Type D at late times [see for 
example [54, 55]]. As a graphic depiction of this evolu- 
tion of spacetime, one may generate PND diagrams sim- 
ilar to Fig. 17 along timclikc worldlincs (instead of null 
geodesies). If one desires a quantitative estimate of the 
rate at which this "settling down" occurs, a metric on the 
anti-celestial sphere has to be defined in order to calcu- 
late distance between the PNDs. However a unique pre- 
scription of this metric requires a unique prescription of 
the tetrad, since Lorentz transformations on tetrad result 
in conformal transformations on the anti-celestial sphere. 
The QKT construction is useful in this context, because 
it uniquely prescribes all the tetrad degrees of freedom si- 
multaneously across spacetime, including along timelike 
worldlines. 



VII. CONCLUSION 

As the numerical relativity codes mature, it is becom- 
ing increasingly important to introduce a protocol that 
allows us to extract and compare physics from these codes 
in an unambiguous fashion. In particular, the quantities 
computed should be independent of the gauge or the for- 
mulation used. Ideally, such a protocol should be valid 
in the strong-field and wave zones and meet the physical 
criteria outlined in Sec. III. 

In this paper, we have suggested one such approach. 



Based on the Newman-Penrose formalism, our method 
fully specifies the tetrad degrees of freedom using purely 
geometric considerations, and two of the gauge degrees 
of freedom are also uniquely fixed using the curvature in- 
variants I and J. In particular, our tetrad construction 
makes use of the quasi-Kinnersley frame (QKF) [9-13], 
which is a transverse frame (TF) that contains the Kin- 
ncrsley tetrad in the Kerr limit. 

By exploiting the relationship between QKF and eigen- 
vectors of the matrix representation Q [as in Eq. (2.12)] 
of the Weyl tensor, one can arrive at several insights re- 
garding the physical properties of the QKF: i) its null 
vector I has a spatial projection pointing along the super- 
Poynting vector [Sec. Ill A] and thus along the direction 
of wave propagation, and ii) there is [Sec. Ill H] a close 
relationship [Fig. 3(b)] between the QKF null basis and 
principal null directions (PNDs) that makes the QKF 
naturally suited to measuring how quickly PNDs bunch 
together (as they converge onto I). These features help 
Newman-Penrose scalars extracted using a QKT to fall 
off correctly in accordance with predictions by the peel- 
ing theorem. 

In the QKF, the eigenvalue if? 2 of the complex matrix 
Q corresponding to the eigenvector that gave us QKF 
is a curvature invariant thus independent of the slicing 
in which the calculation was performed. The physical 
interpretation of ^2 is that it represents the Coulomb 
background portion of Weyl tensor; using this quan- 
tity, we define a pair of geometric coordinates f and 9 
[Sec. HID]. These geometric coordinates vividly depict 
the multipolar structure in the Coulomb potential (as 
can be seen from Figs. 10 and 12). For example, they 
were used to demonstrate that far enough away from 
their source [see Sec. VI B for an empirical cutoff], the 
Coulomb background if?2 [Fig. 10] appears to propagate 
with an almost invariant form along preferred character- 
istic surfaces whose generators are geodesies started off in 
the QKT I direction. Besides fixing the gauge freedom, 
we have also used the differentials df and dO [Sec. HIE] 
to eliminate the spin-boost freedom remaining in QKF, 
yielding a final, gauge-invariant quasi-Kinnersley tetrad 
(QKT). 

As our QKT is constructed from the gauge-invariant 
characteristic structure of Weyl tensor, it can be used 
to explore the physical features of numerical spacetimcs 
in a gauge-invariant way. We have demonstrated this 
desirable property of our QKT with i) a stationary black 
hole spacetime where the hole is offset from the origin 
[Sec. VAl] as well as ii) a gauge wave [Sec. VA2] or 
iii) physical wave [Sec. VB] added to the spacetime of a 
Schwarzschild black hole. These examples serve as useful 
test beds for codes seeking to unambiguously extract the 
physically real effects as opposed to gauge induced false 
signals. 

We have also used the QKT to analyze two equal- 
mass, nonspinning binary-black- hole merger simulations. 
In the first, the two black holes inspiral [Sec. VIA] to- 
ward each other, while in the second they plunge head-on 
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[Sec. VI B]. We have confirmed that the Newman- Penrose 
scalars under the QKT do indeed fall off at the rates ex- 
pected from peeling theorem in these simulations [Fig. 11 
and Fig. 15], and we have explicitly examined the special 
peeling behavior along "critical" directions whose exis- 
tence is ensured by topology [Fig. 17]. 

The gauge invariant feature of the proposed frame- 
work lends itself to several uses. One possible applica- 
tion is that they could help eliminate ambiguities such 
as the pole direction of harmonics used to express grav- 
itational waves [see Sec. VI B 2]. A further application 
is to use the QKT to reduce the ambiguity in measur- 
ing how quickly a spacetime settles down to a Type-D 
spacetime [see the end of Sec. VI B]. The QKT also is 
promising as a wave extraction method that can be per- 
formed real time and ensures that waveform approaches 
its asymptotic value at infinity as rapidly as possible (this 
is illustrated in Fig. 15). For future work, we plan to 
make a comparison between QKT-based wave extraction 
and other wave extraction techniques (such as Cauchy 
Characteristic Extraction [56-61]) using various numeri- 
cal simulations with generic initial conditions. 

The validity of geometric coordinates and the QKT 
throughout spacetime, including the strong field regions, 
suggests that they could be effectively utilized as a vi- 
sualization and diagnostic tool capable of tracking the 
evolution of dynamical features of the spacetime. For 
example, in Fig. 4 (a), the geometric coordinates point 
out the missing rotating quadrupolar moment in the ini- 
tial data. We expect this type of visualization to prove 
valuable in ongoing efforts to construct more realistic ini- 
tial data and reduce spurious "junk" radiation. Other 
utilities for the geometric coordinates include, e.g., their 
potential for helping to improve boundary matching algo- 



rithms [see Fig. 4 (b)]. We would further like to examine 
in greater depth, with the help of QKT, the mechanics 
behind the changes in waveform, as one moves closer to 
the source region [see Fig. 15]. 

Finally, we note that there should exist a close rela- 
tionship between the QKT and the tendex and vortex 
infrastructure introduced in [16, 17], which is based on 
the real eigenvectors and eigenvalues of £ and B. Our 
geometric coordinates and QKT are, in contrast, based 
on the complex eigenvalues and eigenvectors of Q. We 
expect this connection to yield important insights such 
as the slicing dependence of the tendexes and vortexes. 
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